Changeset 2385 for trunk/Clp/src/ClpMatrixBase.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (4 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpMatrixBase.cpp
r2235 r2385 20 20 // Default Constructor 21 21 // 22 ClpMatrixBase::ClpMatrixBase () : 23 rhsOffset_(NULL), 24 startFraction_(0.0), 25 endFraction_(1.0), 26 savedBestDj_(0.0), 27 originalWanted_(0), 28 currentWanted_(0), 29 savedBestSequence_(1), 30 type_(1), 31 lastRefresh_(1), 32 refreshFrequency_(0), 33 minimumObjectsScan_(1), 34 minimumGoodReducedCosts_(1), 35 trueSequenceIn_(1), 36 trueSequenceOut_(1), 37 skipDualCheck_(false) 38 { 39 22 ClpMatrixBase::ClpMatrixBase() 23 : rhsOffset_(NULL) 24 , startFraction_(0.0) 25 , endFraction_(1.0) 26 , savedBestDj_(0.0) 27 , originalWanted_(0) 28 , currentWanted_(0) 29 , savedBestSequence_(1) 30 , type_(1) 31 , lastRefresh_(1) 32 , refreshFrequency_(0) 33 , minimumObjectsScan_(1) 34 , minimumGoodReducedCosts_(1) 35 , trueSequenceIn_(1) 36 , trueSequenceOut_(1) 37 , skipDualCheck_(false) 38 { 40 39 } 41 40 … … 43 42 // Copy constructor 44 43 // 45 ClpMatrixBase::ClpMatrixBase (const ClpMatrixBase & rhs) :46 type_(rhs.type_),47 48 { 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 44 ClpMatrixBase::ClpMatrixBase(const ClpMatrixBase &rhs) 45 : type_(rhs.type_) 46 , skipDualCheck_(rhs.skipDualCheck_) 47 { 48 startFraction_ = rhs.startFraction_; 49 endFraction_ = rhs.endFraction_; 50 savedBestDj_ = rhs.savedBestDj_; 51 originalWanted_ = rhs.originalWanted_; 52 currentWanted_ = rhs.currentWanted_; 53 savedBestSequence_ = rhs.savedBestSequence_; 54 lastRefresh_ = rhs.lastRefresh_; 55 refreshFrequency_ = rhs.refreshFrequency_; 56 minimumObjectsScan_ = rhs.minimumObjectsScan_; 57 minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_; 58 trueSequenceIn_ = rhs.trueSequenceIn_; 59 trueSequenceOut_ = rhs.trueSequenceOut_; 60 skipDualCheck_ = rhs.skipDualCheck_; 61 int numberRows = rhs.getNumRows(); 62 if (rhs.rhsOffset_ && numberRows) { 63 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows); 64 } else { 65 rhsOffset_ = NULL; 66 } 68 67 } 69 68 … … 71 70 // Destructor 72 71 // 73 ClpMatrixBase::~ClpMatrixBase 74 { 75 delete[] rhsOffset_;72 ClpMatrixBase::~ClpMatrixBase() 73 { 74 delete[] rhsOffset_; 76 75 } 77 76 … … 80 79 // 81 80 ClpMatrixBase & 82 ClpMatrixBase::operator=(const ClpMatrixBase &rhs)83 { 84 85 86 delete[] rhsOffset_;87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 81 ClpMatrixBase::operator=(const ClpMatrixBase &rhs) 82 { 83 if (this != &rhs) { 84 type_ = rhs.type_; 85 delete[] rhsOffset_; 86 int numberRows = rhs.getNumRows(); 87 if (rhs.rhsOffset_ && numberRows) { 88 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows); 89 } else { 90 rhsOffset_ = NULL; 91 } 92 startFraction_ = rhs.startFraction_; 93 endFraction_ = rhs.endFraction_; 94 savedBestDj_ = rhs.savedBestDj_; 95 originalWanted_ = rhs.originalWanted_; 96 currentWanted_ = rhs.currentWanted_; 97 savedBestSequence_ = rhs.savedBestSequence_; 98 lastRefresh_ = rhs.lastRefresh_; 99 refreshFrequency_ = rhs.refreshFrequency_; 100 minimumObjectsScan_ = rhs.minimumObjectsScan_; 101 minimumGoodReducedCosts_ = rhs.minimumGoodReducedCosts_; 102 trueSequenceIn_ = rhs.trueSequenceIn_; 103 trueSequenceOut_ = rhs.trueSequenceOut_; 104 skipDualCheck_ = rhs.skipDualCheck_; 105 } 106 return *this; 108 107 } 109 108 // And for scaling  default aborts for when scaling not supported 110 void 111 ClpMatrixBase::times(double scalar, 112 const double * x, double * y, 113 const double * rowScale, 114 const double * /*columnScale*/) const 115 { 116 if (rowScale) { 117 std::cerr << "Scaling not supported  ClpMatrixBase" << std::endl; 118 abort(); 119 } else { 120 times(scalar, x, y); 121 } 109 void ClpMatrixBase::times(double scalar, 110 const double *x, double *y, 111 const double *rowScale, 112 const double * /*columnScale*/) const 113 { 114 if (rowScale) { 115 std::cerr << "Scaling not supported  ClpMatrixBase" << std::endl; 116 abort(); 117 } else { 118 times(scalar, x, y); 119 } 122 120 } 123 121 // And for scaling  default aborts for when scaling not supported 124 void 125 ClpMatrixBase::transposeTimes(double scalar, 126 const double * x, double * y, 127 const double * rowScale, 128 const double * /*columnScale*/, 129 double * /*spare*/) const 130 { 131 if (rowScale) { 132 std::cerr << "Scaling not supported  ClpMatrixBase" << std::endl; 133 abort(); 134 } else { 135 transposeTimes(scalar, x, y); 136 } 122 void ClpMatrixBase::transposeTimes(double scalar, 123 const double *x, double *y, 124 const double *rowScale, 125 const double * /*columnScale*/, 126 double * /*spare*/) const 127 { 128 if (rowScale) { 129 std::cerr << "Scaling not supported  ClpMatrixBase" << std::endl; 130 abort(); 131 } else { 132 transposeTimes(scalar, x, y); 133 } 137 134 } 138 135 /* Subset clone (without gaps). Duplicates are allowed … … 141 138 sense */ 142 139 ClpMatrixBase * 143 ClpMatrixBase::subsetClone ( 144 int /*numberRows*/, const int * /*whichRows*/, 145 int /*numberColumns*/, const int * /*whichColumns*/) const 146 147 148 { 149 std::cerr << "subsetClone not supported  ClpMatrixBase" << std::endl; 150 abort(); 151 return NULL; 140 ClpMatrixBase::subsetClone( 141 int /*numberRows*/, const int * /*whichRows*/, 142 int /*numberColumns*/, const int * /*whichColumns*/) const 143 144 { 145 std::cerr << "subsetClone not supported  ClpMatrixBase" << std::endl; 146 abort(); 147 return NULL; 152 148 } 153 149 /* Given positive integer weights for each row fills in sum of weights … … 157 153 */ 158 154 CoinBigIndex * 159 ClpMatrixBase::dubiousWeights(const ClpSimplex * 160 { 161 162 CoinBigIndex *weights = new CoinBigIndex[number];163 164 165 166 155 ClpMatrixBase::dubiousWeights(const ClpSimplex *model, int * /*inputWeights*/) const 156 { 157 int number = model>numberRows() + model>numberColumns(); 158 CoinBigIndex *weights = new CoinBigIndex[number]; 159 int i; 160 for (i = 0; i < number; i++) 161 weights[i] = 1; 162 return weights; 167 163 } 168 164 #ifndef CLP_NO_VECTOR 169 165 // Append Columns 170 void 171 ClpMatrixBase::appendCols(int /*number*/, 172 const CoinPackedVectorBase * const * /*columns*/) 173 { 174 std::cerr << "appendCols not supported  ClpMatrixBase" << std::endl; 175 abort(); 166 void ClpMatrixBase::appendCols(int /*number*/, 167 const CoinPackedVectorBase *const * /*columns*/) 168 { 169 std::cerr << "appendCols not supported  ClpMatrixBase" << std::endl; 170 abort(); 176 171 } 177 172 // Append Rows 178 void 179 ClpMatrixBase::appendRows(int /*number*/, 180 const CoinPackedVectorBase * const * /*rows*/) 181 { 182 std::cerr << "appendRows not supported  ClpMatrixBase" << std::endl; 183 abort(); 173 void ClpMatrixBase::appendRows(int /*number*/, 174 const CoinPackedVectorBase *const * /*rows*/) 175 { 176 std::cerr << "appendRows not supported  ClpMatrixBase" << std::endl; 177 abort(); 184 178 } 185 179 #endif … … 187 181 Largest refers to largest absolute value. 188 182 */ 189 void 190 ClpMatrixBase::rangeOfElements(double & smallestNegative, double & largestNegative, 191 double & smallestPositive, double & largestPositive) 192 { 193 smallestNegative = 0.0; 194 largestNegative = 0.0; 195 smallestPositive = 0.0; 196 largestPositive = 0.0; 183 void ClpMatrixBase::rangeOfElements(double &smallestNegative, double &largestNegative, 184 double &smallestPositive, double &largestPositive) 185 { 186 smallestNegative = 0.0; 187 largestNegative = 0.0; 188 smallestPositive = 0.0; 189 largestPositive = 0.0; 197 190 } 198 191 /* The length of a majordimension vector. */ 199 int 200 ClpMatrixBase::getVectorLength(int index) const 201 { 202 return getVectorLengths()[index]; 192 int ClpMatrixBase::getVectorLength(int index) const 193 { 194 return getVectorLengths()[index]; 203 195 } 204 196 // Says whether it can do partial pricing 205 bool 206 ClpMatrixBase::canDoPartialPricing() const 207 { 208 return false; // default is no 197 bool ClpMatrixBase::canDoPartialPricing() const 198 { 199 return false; // default is no 209 200 } 210 201 /* Return <code>x *A</code> in <code>z</code> but … … 212 203 Default cheats with fake CoinIndexedVector and 213 204 then calls subsetTransposeTimes */ 214 void 215 ClpMatrixBase::listTransposeTimes(const ClpSimplex * model, 216 double * x, 217 int * y, 218 int number, 219 double * z) const 220 { 221 CoinIndexedVector pi; 222 CoinIndexedVector list; 223 CoinIndexedVector output; 224 int * saveIndices = list.getIndices(); 225 list.setNumElements(number); 226 list.setIndexVector(y); 227 double * savePi = pi.denseVector(); 228 pi.setDenseVector(x); 229 double * saveOutput = output.denseVector(); 230 output.setDenseVector(z); 231 output.setPacked(); 232 subsetTransposeTimes(model, &pi, &list, &output); 233 // restore settings 234 list.setIndexVector(saveIndices); 235 pi.setDenseVector(savePi); 236 output.setDenseVector(saveOutput); 205 void ClpMatrixBase::listTransposeTimes(const ClpSimplex *model, 206 double *x, 207 int *y, 208 int number, 209 double *z) const 210 { 211 CoinIndexedVector pi; 212 CoinIndexedVector list; 213 CoinIndexedVector output; 214 int *saveIndices = list.getIndices(); 215 list.setNumElements(number); 216 list.setIndexVector(y); 217 double *savePi = pi.denseVector(); 218 pi.setDenseVector(x); 219 double *saveOutput = output.denseVector(); 220 output.setDenseVector(z); 221 output.setPacked(); 222 subsetTransposeTimes(model, &pi, &list, &output); 223 // restore settings 224 list.setIndexVector(saveIndices); 225 pi.setDenseVector(savePi); 226 output.setDenseVector(saveOutput); 237 227 } 238 228 // Partial pricing 239 void 240 ClpMatrixBase::partialPricing(ClpSimplex * , double , double , 241 int & , int & ) 242 { 243 std::cerr << "partialPricing not supported  ClpMatrixBase" << std::endl; 244 abort(); 229 void ClpMatrixBase::partialPricing(ClpSimplex *, double, double, 230 int &, int &) 231 { 232 std::cerr << "partialPricing not supported  ClpMatrixBase" << std::endl; 233 abort(); 245 234 } 246 235 /* expands an updated column to allow for extra rows which the main … … 250 239 This will normally be a noop  it is in for GUB! 251 240 */ 252 int 253 ClpMatrixBase::extendUpdated(ClpSimplex * , CoinIndexedVector * , int ) 254 { 255 return 0; 241 int ClpMatrixBase::extendUpdated(ClpSimplex *, CoinIndexedVector *, int) 242 { 243 return 0; 256 244 } 257 245 /* … … 260 248 Remember to update here when settled down 261 249 */ 262 void 263 ClpMatrixBase::primalExpanded(ClpSimplex * , int ) 250 void ClpMatrixBase::primalExpanded(ClpSimplex *, int) 264 251 { 265 252 } … … 269 256 Remember to update here when settled down 270 257 */ 271 void 272 ClpMatrixBase::dualExpanded(ClpSimplex * , 273 CoinIndexedVector * , 274 double * , int ) 258 void ClpMatrixBase::dualExpanded(ClpSimplex *, 259 CoinIndexedVector *, 260 double *, int) 275 261 { 276 262 } … … 280 266 Remember to update here when settled down 281 267 */ 282 int 283 ClpMatrixBase::generalExpanded(ClpSimplex * model, int mode, int &number) 284 { 285 int returnCode = 0; 286 switch (mode) { 287 // Fill in pivotVariable but not for key variables 288 case 0: { 289 int i; 290 int numberBasic = number; 291 int numberColumns = model>numberColumns(); 292 // Use different array so can build from true pivotVariable_ 293 //int * pivotVariable = model>pivotVariable(); 294 int * pivotVariable = model>rowArray(0)>getIndices(); 295 for (i = 0; i < numberColumns; i++) { 296 if (model>getColumnStatus(i) == ClpSimplex::basic) 297 pivotVariable[numberBasic++] = i; 298 } 299 number = numberBasic; 300 } 301 break; 302 // Do initial extra rows + maximum basic 303 case 2: { 304 number = model>numberRows(); 305 } 306 break; 307 // To see if can dual or primal 308 case 4: { 309 returnCode = 3; 310 } 311 break; 312 default: 313 break; 314 } 315 return returnCode; 268 int ClpMatrixBase::generalExpanded(ClpSimplex *model, int mode, int &number) 269 { 270 int returnCode = 0; 271 switch (mode) { 272 // Fill in pivotVariable but not for key variables 273 case 0: { 274 int i; 275 int numberBasic = number; 276 int numberColumns = model>numberColumns(); 277 // Use different array so can build from true pivotVariable_ 278 //int * pivotVariable = model>pivotVariable(); 279 int *pivotVariable = model>rowArray(0)>getIndices(); 280 for (i = 0; i < numberColumns; i++) { 281 if (model>getColumnStatus(i) == ClpSimplex::basic) 282 pivotVariable[numberBasic++] = i; 283 } 284 number = numberBasic; 285 } break; 286 // Do initial extra rows + maximum basic 287 case 2: { 288 number = model>numberRows(); 289 } break; 290 // To see if can dual or primal 291 case 4: { 292 returnCode = 3; 293 } break; 294 default: 295 break; 296 } 297 return returnCode; 316 298 } 317 299 // Sets up an effective RHS 318 void 319 ClpMatrixBase::useEffectiveRhs(ClpSimplex * ) 320 { 321 std::cerr << "useEffectiveRhs not supported  ClpMatrixBase" << std::endl; 322 abort(); 300 void ClpMatrixBase::useEffectiveRhs(ClpSimplex *) 301 { 302 std::cerr << "useEffectiveRhs not supported  ClpMatrixBase" << std::endl; 303 abort(); 323 304 } 324 305 /* Returns effective RHS if it is being used. This is used for long problems … … 326 307 expensive. This may recompute */ 327 308 double * 328 ClpMatrixBase::rhsOffset(ClpSimplex * 309 ClpMatrixBase::rhsOffset(ClpSimplex *model, bool forceRefresh, bool 329 310 #ifdef CLP_DEBUG 330 check311 check 331 312 #endif 332 333 { 334 313 ) 314 { 315 if (rhsOffset_) { 335 316 #ifdef CLP_DEBUG 336 337 338 339 340 341 double * solution = new double[numberColumns];342 double *rhs = new double[numberRows];343 const double *solutionSlack = model>solutionRegion(0);344 345 346 347 348 349 350 351 352 353 354 355 356 357 delete[] solution;358 359 360 361 362 317 if (check) { 318 // no need  but check anyway 319 // zero out basic 320 int numberRows = model>numberRows(); 321 int numberColumns = model>numberColumns(); 322 double *solution = new double[numberColumns]; 323 double *rhs = new double[numberRows]; 324 const double *solutionSlack = model>solutionRegion(0); 325 CoinMemcpyN(model>solutionRegion(), numberColumns, solution); 326 int iRow; 327 for (iRow = 0; iRow < numberRows; iRow++) { 328 if (model>getRowStatus(iRow) != ClpSimplex::basic) 329 rhs[iRow] = solutionSlack[iRow]; 330 else 331 rhs[iRow] = 0.0; 332 } 333 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 334 if (model>getColumnStatus(iColumn) == ClpSimplex::basic) 335 solution[iColumn] = 0.0; 336 } 337 times(1.0, solution, rhs); 338 delete[] solution; 339 for (iRow = 0; iRow < numberRows; iRow++) { 340 if (fabs(rhs[iRow]  rhsOffset_[iRow]) > 1.0e3) 341 printf("** bad effective %d  true %g old %g\n", iRow, rhs[iRow], rhsOffset_[iRow]); 342 } 343 } 363 344 #endif 364 if (forceRefresh  (refreshFrequency_ && model>numberIterations() >= 365 lastRefresh_ + refreshFrequency_)) { 366 // zero out basic 367 int numberRows = model>numberRows(); 368 int numberColumns = model>numberColumns(); 369 double * solution = new double [numberColumns]; 370 const double * solutionSlack = model>solutionRegion(0); 371 CoinMemcpyN(model>solutionRegion(), numberColumns, solution); 372 for (int iRow = 0; iRow < numberRows; iRow++) { 373 if (model>getRowStatus(iRow) != ClpSimplex::basic) 374 rhsOffset_[iRow] = solutionSlack[iRow]; 375 else 376 rhsOffset_[iRow] = 0.0; 377 } 378 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 379 if (model>getColumnStatus(iColumn) == ClpSimplex::basic) 380 solution[iColumn] = 0.0; 381 } 382 times(1.0, solution, rhsOffset_); 383 delete [] solution; 384 lastRefresh_ = model>numberIterations(); 385 } 386 } 387 return rhsOffset_; 345 if (forceRefresh  (refreshFrequency_ && model>numberIterations() >= lastRefresh_ + refreshFrequency_)) { 346 // zero out basic 347 int numberRows = model>numberRows(); 348 int numberColumns = model>numberColumns(); 349 double *solution = new double[numberColumns]; 350 const double *solutionSlack = model>solutionRegion(0); 351 CoinMemcpyN(model>solutionRegion(), numberColumns, solution); 352 for (int iRow = 0; iRow < numberRows; iRow++) { 353 if (model>getRowStatus(iRow) != ClpSimplex::basic) 354 rhsOffset_[iRow] = solutionSlack[iRow]; 355 else 356 rhsOffset_[iRow] = 0.0; 357 } 358 for (int iColumn = 0; iColumn < numberColumns; iColumn++) { 359 if (model>getColumnStatus(iColumn) == ClpSimplex::basic) 360 solution[iColumn] = 0.0; 361 } 362 times(1.0, solution, rhsOffset_); 363 delete[] solution; 364 lastRefresh_ = model>numberIterations(); 365 } 366 } 367 return rhsOffset_; 388 368 } 389 369 /* 390 370 update information for a pivot (and effective rhs) 391 371 */ 392 int 393 ClpMatrixBase::updatePivot(ClpSimplex * model, double oldInValue, double ) 394 { 395 if (rhsOffset_) { 396 // update effective rhs 397 int sequenceIn = model>sequenceIn(); 398 int sequenceOut = model>sequenceOut(); 399 double * solution = model>solutionRegion(); 400 int numberColumns = model>numberColumns(); 401 if (sequenceIn == sequenceOut) { 402 if (sequenceIn < numberColumns) 403 add(model, rhsOffset_, sequenceIn, oldInValue  solution[sequenceIn]); 404 } else { 405 if (sequenceIn < numberColumns) 406 add(model, rhsOffset_, sequenceIn, oldInValue); 407 if (sequenceOut < numberColumns) 408 add(model, rhsOffset_, sequenceOut, solution[sequenceOut]); 409 } 410 } 411 return 0; 412 } 413 int 414 ClpMatrixBase::hiddenRows() const 415 { 416 return 0; 372 int ClpMatrixBase::updatePivot(ClpSimplex *model, double oldInValue, double) 373 { 374 if (rhsOffset_) { 375 // update effective rhs 376 int sequenceIn = model>sequenceIn(); 377 int sequenceOut = model>sequenceOut(); 378 double *solution = model>solutionRegion(); 379 int numberColumns = model>numberColumns(); 380 if (sequenceIn == sequenceOut) { 381 if (sequenceIn < numberColumns) 382 add(model, rhsOffset_, sequenceIn, oldInValue  solution[sequenceIn]); 383 } else { 384 if (sequenceIn < numberColumns) 385 add(model, rhsOffset_, sequenceIn, oldInValue); 386 if (sequenceOut < numberColumns) 387 add(model, rhsOffset_, sequenceOut, solution[sequenceOut]); 388 } 389 } 390 return 0; 391 } 392 int ClpMatrixBase::hiddenRows() const 393 { 394 return 0; 417 395 } 418 396 /* Creates a variable. This is called after partial pricing and may modify matrix. 419 397 May update bestSequence. 420 398 */ 421 void 422 ClpMatrixBase::createVariable(ClpSimplex *, int &) 399 void ClpMatrixBase::createVariable(ClpSimplex *, int &) 423 400 { 424 401 } 425 402 // Returns reduced cost of a variable 426 403 double 427 ClpMatrixBase::reducedCost(ClpSimplex * 428 { 429 430 431 432 433 434 404 ClpMatrixBase::reducedCost(ClpSimplex *model, int sequence) const 405 { 406 int numberRows = model>numberRows(); 407 int numberColumns = model>numberColumns(); 408 if (sequence < numberRows + numberColumns) 409 return model>djRegion()[sequence]; 410 else 411 return savedBestDj_; 435 412 } 436 413 /* Just for debug if odd type matrix. 437 414 Returns number and sum of primal infeasibilities. 438 415 */ 439 int 440 ClpMatrixBase::checkFeasible(ClpSimplex * model, double & sum) const 441 { 442 int numberRows = model>numberRows(); 443 double * rhs = new double[numberRows]; 444 int numberColumns = model>numberColumns(); 445 int iRow; 446 CoinZeroN(rhs, numberRows); 447 times(1.0, model>solutionRegion(), rhs, model>rowScale(), model>columnScale()); 448 int iColumn; 449 int logLevel = model>messageHandler()>logLevel(); 450 int numberInfeasible = 0; 451 const double * rowLower = model>lowerRegion(0); 452 const double * rowUpper = model>upperRegion(0); 453 const double * solution; 454 solution = model>solutionRegion(0); 455 double tolerance = model>primalTolerance() * 1.01; 456 sum = 0.0; 457 for (iRow = 0; iRow < numberRows; iRow++) { 458 double value = rhs[iRow]; 459 double value2 = solution[iRow]; 460 if (logLevel > 3) { 461 if (fabs(value  value2) > 1.0e8) 462 printf("Row %d stored %g, computed %g\n", iRow, value2, value); 463 } 464 if (value < rowLower[iRow]  tolerance  465 value > rowUpper[iRow] + tolerance) { 466 numberInfeasible++; 467 sum += CoinMax(rowLower[iRow]  value, value  rowUpper[iRow]); 468 } 469 if (value2 > rowLower[iRow] + tolerance && 470 value2 < rowUpper[iRow]  tolerance && 471 model>getRowStatus(iRow) != ClpSimplex::basic) { 472 assert (model>getRowStatus(iRow) == ClpSimplex::superBasic); 473 } 474 } 475 const double * columnLower = model>lowerRegion(1); 476 const double * columnUpper = model>upperRegion(1); 477 solution = model>solutionRegion(1); 478 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 479 double value = solution[iColumn]; 480 if (value < columnLower[iColumn]  tolerance  481 value > columnUpper[iColumn] + tolerance) { 482 numberInfeasible++; 483 sum += CoinMax(columnLower[iColumn]  value, value  columnUpper[iColumn]); 484 } 485 if (value > columnLower[iColumn] + tolerance && 486 value < columnUpper[iColumn]  tolerance && 487 model>getColumnStatus(iColumn) != ClpSimplex::basic) { 488 assert (model>getColumnStatus(iColumn) == ClpSimplex::superBasic); 489 } 490 } 491 delete [] rhs; 492 return numberInfeasible; 416 int ClpMatrixBase::checkFeasible(ClpSimplex *model, double &sum) const 417 { 418 int numberRows = model>numberRows(); 419 double *rhs = new double[numberRows]; 420 int numberColumns = model>numberColumns(); 421 int iRow; 422 CoinZeroN(rhs, numberRows); 423 times(1.0, model>solutionRegion(), rhs, model>rowScale(), model>columnScale()); 424 int iColumn; 425 int logLevel = model>messageHandler()>logLevel(); 426 int numberInfeasible = 0; 427 const double *rowLower = model>lowerRegion(0); 428 const double *rowUpper = model>upperRegion(0); 429 const double *solution; 430 solution = model>solutionRegion(0); 431 double tolerance = model>primalTolerance() * 1.01; 432 sum = 0.0; 433 for (iRow = 0; iRow < numberRows; iRow++) { 434 double value = rhs[iRow]; 435 double value2 = solution[iRow]; 436 if (logLevel > 3) { 437 if (fabs(value  value2) > 1.0e8) 438 printf("Row %d stored %g, computed %g\n", iRow, value2, value); 439 } 440 if (value < rowLower[iRow]  tolerance  value > rowUpper[iRow] + tolerance) { 441 numberInfeasible++; 442 sum += CoinMax(rowLower[iRow]  value, value  rowUpper[iRow]); 443 } 444 if (value2 > rowLower[iRow] + tolerance && value2 < rowUpper[iRow]  tolerance && model>getRowStatus(iRow) != ClpSimplex::basic) { 445 assert(model>getRowStatus(iRow) == ClpSimplex::superBasic); 446 } 447 } 448 const double *columnLower = model>lowerRegion(1); 449 const double *columnUpper = model>upperRegion(1); 450 solution = model>solutionRegion(1); 451 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 452 double value = solution[iColumn]; 453 if (value < columnLower[iColumn]  tolerance  value > columnUpper[iColumn] + tolerance) { 454 numberInfeasible++; 455 sum += CoinMax(columnLower[iColumn]  value, value  columnUpper[iColumn]); 456 } 457 if (value > columnLower[iColumn] + tolerance && value < columnUpper[iColumn]  tolerance && model>getColumnStatus(iColumn) != ClpSimplex::basic) { 458 assert(model>getColumnStatus(iColumn) == ClpSimplex::superBasic); 459 } 460 } 461 delete[] rhs; 462 return numberInfeasible; 493 463 } 494 464 // These have to match ClpPrimalColumnSteepest version 495 #define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)465 #define reference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0) 496 466 // Updates second array for steepest and does devex weights (need not be coded) 497 void 498 ClpMatrixBase::subsetTimes2(const ClpSimplex * model, 499 CoinIndexedVector * dj1, 500 const CoinIndexedVector * pi2, CoinIndexedVector * dj2, 501 double referenceIn, double devex, 502 // Array for exact devex to say what is in reference framework 503 unsigned int * reference, 504 double * weights, double scaleFactor) 505 { 506 // get subset which have nonzero tableau elements 507 subsetTransposeTimes(model, pi2, dj1, dj2); 508 bool killDjs = (scaleFactor == 0.0); 509 if (!scaleFactor) 510 scaleFactor = 1.0; 511 // columns 512 513 int number = dj1>getNumElements(); 514 const int * index = dj1>getIndices(); 515 double * updateBy = dj1>denseVector(); 516 double * updateBy2 = dj2>denseVector(); 517 518 for (int j = 0; j < number; j++) { 519 double thisWeight; 520 double pivot; 521 double pivotSquared; 522 int iSequence = index[j]; 523 double value2 = updateBy[j]; 524 if (killDjs) 525 updateBy[j] = 0.0; 526 double modification = updateBy2[j]; 527 updateBy2[j] = 0.0; 528 ClpSimplex::Status status = model>getStatus(iSequence); 529 530 if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) { 531 thisWeight = weights[iSequence]; 532 pivot = value2 * scaleFactor; 533 pivotSquared = pivot * pivot; 534 535 thisWeight += pivotSquared * devex + pivot * modification; 536 if (thisWeight < DEVEX_TRY_NORM) { 537 if (referenceIn < 0.0) { 538 // steepest 539 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 540 } else { 541 // exact 542 thisWeight = referenceIn * pivotSquared; 543 if (reference(iSequence)) 544 thisWeight += 1.0; 545 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 546 } 547 } 548 weights[iSequence] = thisWeight; 549 } 550 } 551 dj2>setNumElements(0); 467 void ClpMatrixBase::subsetTimes2(const ClpSimplex *model, 468 CoinIndexedVector *dj1, 469 const CoinIndexedVector *pi2, CoinIndexedVector *dj2, 470 double referenceIn, double devex, 471 // Array for exact devex to say what is in reference framework 472 unsigned int *reference, 473 double *weights, double scaleFactor) 474 { 475 // get subset which have nonzero tableau elements 476 subsetTransposeTimes(model, pi2, dj1, dj2); 477 bool killDjs = (scaleFactor == 0.0); 478 if (!scaleFactor) 479 scaleFactor = 1.0; 480 // columns 481 482 int number = dj1>getNumElements(); 483 const int *index = dj1>getIndices(); 484 double *updateBy = dj1>denseVector(); 485 double *updateBy2 = dj2>denseVector(); 486 487 for (int j = 0; j < number; j++) { 488 double thisWeight; 489 double pivot; 490 double pivotSquared; 491 int iSequence = index[j]; 492 double value2 = updateBy[j]; 493 if (killDjs) 494 updateBy[j] = 0.0; 495 double modification = updateBy2[j]; 496 updateBy2[j] = 0.0; 497 ClpSimplex::Status status = model>getStatus(iSequence); 498 499 if (status != ClpSimplex::basic && status != ClpSimplex::isFixed) { 500 thisWeight = weights[iSequence]; 501 pivot = value2 * scaleFactor; 502 pivotSquared = pivot * pivot; 503 504 thisWeight += pivotSquared * devex + pivot * modification; 505 if (thisWeight < DEVEX_TRY_NORM) { 506 if (referenceIn < 0.0) { 507 // steepest 508 thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); 509 } else { 510 // exact 511 thisWeight = referenceIn * pivotSquared; 512 if (reference(iSequence)) 513 thisWeight += 1.0; 514 thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); 515 } 516 } 517 weights[iSequence] = thisWeight; 518 } 519 } 520 dj2>setNumElements(0); 552 521 } 553 522 // Correct sequence in and out to give true value 554 void 555 ClpMatrixBase::correctSequence(const ClpSimplex * , int & , int & ) 523 void ClpMatrixBase::correctSequence(const ClpSimplex *, int &, int &) 556 524 { 557 525 } 558 526 // Really scale matrix 559 void 560 ClpMatrixBase::reallyScale(const double * , const double * ) 561 { 562 std::cerr << "reallyScale not supported  ClpMatrixBase" << std::endl; 563 abort(); 527 void ClpMatrixBase::reallyScale(const double *, const double *) 528 { 529 std::cerr << "reallyScale not supported  ClpMatrixBase" << std::endl; 530 abort(); 564 531 } 565 532 // Updates two arrays for steepest 566 int 567 ClpMatrixBase::transposeTimes2(const ClpSimplex * , 568 const CoinIndexedVector * , CoinIndexedVector *, 569 const CoinIndexedVector * , 570 CoinIndexedVector * , 571 double * , double * , 572 double , double , 573 // Array for exact devex to say what is in reference framework 574 unsigned int * , 575 double * , double ) 576 { 577 std::cerr << "transposeTimes2 not supported  ClpMatrixBase" << std::endl; 578 abort(); 579 return 0; 533 int ClpMatrixBase::transposeTimes2(const ClpSimplex *, 534 const CoinIndexedVector *, CoinIndexedVector *, 535 const CoinIndexedVector *, 536 CoinIndexedVector *, 537 double *, double *, 538 double, double, 539 // Array for exact devex to say what is in reference framework 540 unsigned int *, 541 double *, double) 542 { 543 std::cerr << "transposeTimes2 not supported  ClpMatrixBase" << std::endl; 544 abort(); 545 return 0; 580 546 } 581 547 /* Set the dimensions of the matrix. In effect, append new empty … … 584 550 MUST be at least as large as the current ones otherwise an exception 585 551 is thrown. */ 586 void 587 ClpMatrixBase::setDimensions(int , int ) 588 { 589 // If odd matrix assume user knows what they are doing 552 void ClpMatrixBase::setDimensions(int, int) 553 { 554 // If odd matrix assume user knows what they are doing 590 555 } 591 556 /* Append a set of rows/columns to the end of the matrix. Returns number of errors … … 593 558 number of columns1/rows1 (if numberOther>0) or duplicates 594 559 If 0 then rows, 1 if columns */ 595 int 596 ClpMatrixBase::appendMatrix(int , int , 597 const CoinBigIndex * , const int * , 598 const double * , int ) 599 { 600 std::cerr << "appendMatrix not supported  ClpMatrixBase" << std::endl; 601 abort(); 602 return 1; 560 int ClpMatrixBase::appendMatrix(int, int, 561 const CoinBigIndex *, const int *, 562 const double *, int) 563 { 564 std::cerr << "appendMatrix not supported  ClpMatrixBase" << std::endl; 565 abort(); 566 return 1; 603 567 } 604 568 … … 606 570 This works for either ordering If the new element is zero it will be 607 571 deleted unless keepZero true */ 608 void 609 ClpMatrixBase::modifyCoefficient(int , int , double , 610 bool ) 611 { 612 std::cerr << "modifyCoefficient not supported  ClpMatrixBase" << std::endl; 613 abort(); 572 void ClpMatrixBase::modifyCoefficient(int, int, double, 573 bool) 574 { 575 std::cerr << "modifyCoefficient not supported  ClpMatrixBase" << std::endl; 576 abort(); 614 577 } 615 578 #if COIN_LONG_WORK 616 579 // For long double versions (aborts if not supported) 617 void 618 ClpMatrixBase::times(CoinWorkDouble scalar, 619 const CoinWorkDouble * x, CoinWorkDouble * y) const 620 { 621 std::cerr << "long times not supported  ClpMatrixBase" << std::endl; 622 abort(); 623 } 624 void 625 ClpMatrixBase::transposeTimes(CoinWorkDouble scalar, 626 const CoinWorkDouble * x, CoinWorkDouble * y) const 627 { 628 std::cerr << "long transposeTimes not supported  ClpMatrixBase" << std::endl; 629 abort(); 580 void ClpMatrixBase::times(CoinWorkDouble scalar, 581 const CoinWorkDouble *x, CoinWorkDouble *y) const 582 { 583 std::cerr << "long times not supported  ClpMatrixBase" << std::endl; 584 abort(); 585 } 586 void ClpMatrixBase::transposeTimes(CoinWorkDouble scalar, 587 const CoinWorkDouble *x, CoinWorkDouble *y) const 588 { 589 std::cerr << "long transposeTimes not supported  ClpMatrixBase" << std::endl; 590 abort(); 630 591 } 631 592 #endif 593 594 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 595 */
Note: See TracChangeset
for help on using the changeset viewer.