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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/AbcSimplexFactorization.cpp
r2166 r2385 52 52 // Default Constructor 53 53 // 54 AbcSimplexFactorization::AbcSimplexFactorization 55 { 56 model_ =NULL;54 AbcSimplexFactorization::AbcSimplexFactorization(int /*numberRows*/) 55 { 56 model_ = NULL; 57 57 coinAbcFactorization_ = new CoinAbcFactorization(); 58 58 forceB_ = 0; … … 60 60 goSmallThreshold_ = USE_SMALL_FAC; 61 61 goLongThreshold_ = USE_LONG_FAC; 62 numberSlacks_ =0;62 numberSlacks_ = 0; 63 63 } 64 64 … … 66 66 // Copy constructor 67 67 // 68 AbcSimplexFactorization::AbcSimplexFactorization (const AbcSimplexFactorization &rhs,69 68 AbcSimplexFactorization::AbcSimplexFactorization(const AbcSimplexFactorization &rhs, 69 int denseIfSmaller) 70 70 { 71 71 forceB_ = rhs.forceB_; … … 73 73 goSmallThreshold_ = rhs.goSmallThreshold_; 74 74 goLongThreshold_ = rhs.goLongThreshold_; 75 numberSlacks_ =rhs.numberSlacks_;76 model_ =rhs.model_;75 numberSlacks_ = rhs.numberSlacks_; 76 model_ = rhs.model_; 77 77 #ifndef ABC_USE_COIN_FACTORIZATION 78 78 int goDense = 0; 79 79 if (denseIfSmaller > 0 && denseIfSmaller <= goDenseThreshold_) { 80 CoinAbcDenseFactorization * denseR = 81 dynamic_cast<CoinAbcDenseFactorization *>(rhs.coinAbcFactorization_); 80 CoinAbcDenseFactorization *denseR = dynamic_cast< CoinAbcDenseFactorization * >(rhs.coinAbcFactorization_); 82 81 if (!denseR) 83 82 goDense = 1; … … 112 111 else 113 112 coinAbcFactorization_ = new CoinAbcFactorization(); 114 assert 113 assert(coinAbcFactorization_); 115 114 coinAbcFactorization_>maximumPivots(rhs.coinAbcFactorization_>maximumPivots()); 116 115 coinAbcFactorization_>pivotTolerance(rhs.coinAbcFactorization_>pivotTolerance()); … … 124 123 // Destructor 125 124 // 126 AbcSimplexFactorization::~AbcSimplexFactorization 125 AbcSimplexFactorization::~AbcSimplexFactorization() 127 126 { 128 127 delete coinAbcFactorization_; … … 133 132 // 134 133 AbcSimplexFactorization & 135 AbcSimplexFactorization::operator=(const AbcSimplexFactorization &rhs)134 AbcSimplexFactorization::operator=(const AbcSimplexFactorization &rhs) 136 135 { 137 136 if (this != &rhs) { 138 137 forceB_ = rhs.forceB_; 139 model_ =rhs.model_;138 model_ = rhs.model_; 140 139 goDenseThreshold_ = rhs.goDenseThreshold_; 141 140 goSmallThreshold_ = rhs.goSmallThreshold_; 142 141 goLongThreshold_ = rhs.goLongThreshold_; 143 numberSlacks_ =rhs.numberSlacks_;144 142 numberSlacks_ = rhs.numberSlacks_; 143 145 144 if (rhs.coinAbcFactorization_) { 146 145 delete coinAbcFactorization_; … … 155 154 } 156 155 // Go over to dense code 157 void 158 AbcSimplexFactorization::goDenseOrSmall(int numberRows) 156 void AbcSimplexFactorization::goDenseOrSmall(int numberRows) 159 157 { 160 158 #ifndef ABC_USE_COIN_FACTORIZATION … … 174 172 } 175 173 // If nonzero force use of 1,dense 2,small 3,long 176 void 177 AbcSimplexFactorization::forceOtherFactorization(int which) 174 void AbcSimplexFactorization::forceOtherFactorization(int which) 178 175 { 179 176 #ifndef ABC_USE_COIN_FACTORIZATION … … 205 202 } 206 203 #ifdef CLP_FACTORIZATION_NEW_TIMING 207 static bool readTwiddle =false;208 static double weightIncU =1.0;209 static double weightR =2.0;210 static double weightRest =1.0;211 static double weightFactL =30.0;212 static double weightFactDense =0.1;213 static double weightNrows =10.0;214 static double increaseNeeded =1.1;204 static bool readTwiddle = false; 205 static double weightIncU = 1.0; 206 static double weightR = 2.0; 207 static double weightRest = 1.0; 208 static double weightFactL = 30.0; 209 static double weightFactDense = 0.1; 210 static double weightNrows = 10.0; 211 static double increaseNeeded = 1.1; 215 212 static double constWeightIterate = 1.0; 216 213 static double weightNrowsIterate = 3.0; 217 bool 218 AbcSimplexFactorization::timeToRefactorize() const 219 { 220 bool reFactor = (coinAbcFactorization_>pivots() * 3 > coinAbcFactorization_>maximumPivots() * 2 && 221 coinAbcFactorization_>numberElementsR() * 3 > (coinAbcFactorization_>numberElementsL() + 222 coinAbcFactorization_>numberElementsU()) * 2 + 1000 && 223 !coinAbcFactorization_>numberDense()); 224 bool reFactor3=false; 225 int numberPivots=coinAbcFactorization_>pivots(); 214 bool AbcSimplexFactorization::timeToRefactorize() const 215 { 216 bool reFactor = (coinAbcFactorization_>pivots() * 3 > coinAbcFactorization_>maximumPivots() * 2 && coinAbcFactorization_>numberElementsR() * 3 > (coinAbcFactorization_>numberElementsL() + coinAbcFactorization_>numberElementsU()) * 2 + 1000 && !coinAbcFactorization_>numberDense()); 217 bool reFactor3 = false; 218 int numberPivots = coinAbcFactorization_>pivots(); 226 219 //if (coinAbcFactorization_>pivots()<2) 227 if (numberPivots >lastNumberPivots_) {220 if (numberPivots > lastNumberPivots_) { 228 221 if (!lastNumberPivots_) { 229 222 //lastR=0; 230 223 //lastU=endLengthU; 231 totalInR_ =0.0;232 totalInIncreasingU_ =0.0;233 shortestAverage_ =COIN_DBL_MAX;224 totalInR_ = 0.0; 225 totalInIncreasingU_ = 0.0; 226 shortestAverage_ = COIN_DBL_MAX; 234 227 if (!readTwiddle) { 235 readTwiddle=true;236 char *environ = getenv("CLP_TWIDDLE");237 238 sscanf(environ,"%lg %lg %lg %lg %lg %lg %lg %lg %lg",239 &weightIncU,&weightR,&weightRest,&weightFactL,240 &weightFactDense,&weightNrows,&increaseNeeded,241 &constWeightIterate,&weightNrowsIterate);242 243 244 weightIncU,weightR,weightRest,weightFactL,245 weightFactDense,weightNrows,increaseNeeded,246 constWeightIterate,weightNrowsIterate);228 readTwiddle = true; 229 char *environ = getenv("CLP_TWIDDLE"); 230 if (environ) { 231 sscanf(environ, "%lg %lg %lg %lg %lg %lg %lg %lg %lg", 232 &weightIncU, &weightR, &weightRest, &weightFactL, 233 &weightFactDense, &weightNrows, &increaseNeeded, 234 &constWeightIterate, &weightNrowsIterate); 235 } 236 printf("weightIncU %g, weightR %g, weightRest %g, weightFactL %g, weightFactDense %g, weightNrows %g increaseNeeded %g constWeightIterate %g weightNrowsIterate %g\n", 237 weightIncU, weightR, weightRest, weightFactL, 238 weightFactDense, weightNrows, increaseNeeded, 239 constWeightIterate, weightNrowsIterate); 247 240 } 248 241 } 249 lastNumberPivots_ =numberPivots;250 int numberDense =coinAbcFactorization_>numberDense();251 double nnd =numberDense*numberDense;252 int lengthL =coinAbcFactorization_>numberElementsL();253 int lengthR =coinAbcFactorization_>numberElementsR();242 lastNumberPivots_ = numberPivots; 243 int numberDense = coinAbcFactorization_>numberDense(); 244 double nnd = numberDense * numberDense; 245 int lengthL = coinAbcFactorization_>numberElementsL(); 246 int lengthR = coinAbcFactorization_>numberElementsR(); 254 247 int numberRows = coinAbcFactorization_>numberRows(); 255 int lengthU=coinAbcFactorization_>numberElementsU() 256 (numberRowsnumberDense); 248 int lengthU = coinAbcFactorization_>numberElementsU()  (numberRows  numberDense); 257 249 totalInR_ += lengthR; 258 int effectiveU =lengthUeffectiveStartNumberU_;250 int effectiveU = lengthU  effectiveStartNumberU_; 259 251 totalInIncreasingU_ += effectiveU; 260 252 //lastR=lengthR; 261 253 //lastU=lengthU; 262 double rest=lengthL+0.05*nnd; 263 double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd 264 + weightNrows*numberRows; 265 double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_) 266 + weightNrowsIterate*numberRows; 267 double variableWeight = weightIncU*totalInIncreasingU_+ 268 weightR*totalInR_+weightRest*rest; 269 double average=constWeightIterateX+ 270 (constWeightFactor+variableWeight)/static_cast<double>(numberPivots); 254 double rest = lengthL + 0.05 * nnd; 255 double constWeightFactor = weightFactL * lengthL + weightFactDense * nnd 256 + weightNrows * numberRows; 257 double constWeightIterateX = constWeightIterate * (lengthL + endLengthU_) 258 + weightNrowsIterate * numberRows; 259 double variableWeight = weightIncU * totalInIncreasingU_ + weightR * totalInR_ + weightRest * rest; 260 double average = constWeightIterateX + (constWeightFactor + variableWeight) / static_cast< double >(numberPivots); 271 261 #if 0 272 262 if ((numberPivots%20)==0&&!ifPrint3) … … 275 265 lengthU,lengthL,lengthR,nnd,average); 276 266 #endif 277 shortestAverage_=CoinMin(shortestAverage_,average); 278 if (average>increaseNeeded*shortestAverage_&& 279 coinAbcFactorization_>pivots()>30) { 267 shortestAverage_ = CoinMin(shortestAverage_, average); 268 if (average > increaseNeeded * shortestAverage_ && coinAbcFactorization_>pivots() > 30) { 280 269 //printf("PIVX %d nrow %d startU %d now %d L %d R %d dense %g average %g\n", 281 270 // numberPivots,numberRows,effectiveStartNumberU_, 282 271 // lengthU,lengthL,lengthR,nnd,average); 283 reFactor3 =true;284 } 285 } 286 if (reFactor  reFactor3) {287 reFactor =true;272 reFactor3 = true; 273 } 274 } 275 if (reFactor  reFactor3) { 276 reFactor = true; 288 277 } 289 278 return reFactor; 290 279 } 291 #if CLP_FACTORIZATION_NEW_TIMING>1 292 void 293 AbcSimplexFactorization::statsRefactor(char when) const 294 { 295 int numberPivots=coinAbcFactorization_>pivots(); 296 int numberDense=coinAbcFactorization_>numberDense(); 297 double nnd=numberDense*numberDense; 298 int lengthL=coinAbcFactorization_>numberElementsL(); 299 int lengthR=coinAbcFactorization_>numberElementsR(); 280 #if CLP_FACTORIZATION_NEW_TIMING > 1 281 void AbcSimplexFactorization::statsRefactor(char when) const 282 { 283 int numberPivots = coinAbcFactorization_>pivots(); 284 int numberDense = coinAbcFactorization_>numberDense(); 285 double nnd = numberDense * numberDense; 286 int lengthL = coinAbcFactorization_>numberElementsL(); 287 int lengthR = coinAbcFactorization_>numberElementsR(); 300 288 int numberRows = coinAbcFactorization_>numberRows(); 301 int lengthU=coinAbcFactorization_>numberElementsU() 302 (numberRowsnumberDense); 303 double rest=lengthL+0.05*nnd; 304 double constWeightFactor = weightFactL*lengthL+weightFactDense*nnd 305 + weightNrows*numberRows; 306 double constWeightIterateX = constWeightIterate*(lengthL+endLengthU_) 307 + weightNrowsIterate*numberRows; 308 double variableWeight = weightIncU*totalInIncreasingU_+ 309 weightR*totalInR_+weightRest*rest; 310 double average=constWeightIterateX+ 311 (constWeightFactor+variableWeight)/static_cast<double>(numberPivots); 289 int lengthU = coinAbcFactorization_>numberElementsU()  (numberRows  numberDense); 290 double rest = lengthL + 0.05 * nnd; 291 double constWeightFactor = weightFactL * lengthL + weightFactDense * nnd 292 + weightNrows * numberRows; 293 double constWeightIterateX = constWeightIterate * (lengthL + endLengthU_) 294 + weightNrowsIterate * numberRows; 295 double variableWeight = weightIncU * totalInIncreasingU_ + weightR * totalInR_ + weightRest * rest; 296 double average = constWeightIterateX + (constWeightFactor + variableWeight) / static_cast< double >(numberPivots); 312 297 printf("APIV%c %d nrow %d startU %d now %d L %d R %d dense %g average %g  shortest %g\n", 313 when,numberPivots,numberRows,effectiveStartNumberU_, 314 lengthU,lengthL,lengthR,nnd,average,shortestAverage_); 315 } 316 #endif 317 #else 318 bool 319 AbcSimplexFactorization::timeToRefactorize() const 298 when, numberPivots, numberRows, effectiveStartNumberU_, 299 lengthU, lengthL, lengthR, nnd, average, shortestAverage_); 300 } 301 #endif 302 #else 303 bool AbcSimplexFactorization::timeToRefactorize() const 320 304 { 321 305 return coinAbcFactorization_>pivots() > coinAbcFactorization_>numberRows() / 2.45 + 20; … … 324 308 /* returns empty fake vector carved out of existing 325 309 later  maybe use associated arrays */ 326 static CoinIndexedVector * fakeVector(CoinIndexedVector *vector,327 328 { 329 int oldCapacity =vector>capacity();330 CoinIndexedVector * 310 static CoinIndexedVector *fakeVector(CoinIndexedVector *vector, 311 int fakeCapacity) 312 { 313 int oldCapacity = vector>capacity(); 314 CoinIndexedVector *newVector = new CoinIndexedVector(); 331 315 newVector>setCapacity(fakeCapacity); 332 newVector>setDenseVector(vector>denseVector()+oldCapacity); 333 newVector>setIndexVector(vector>getIndices()+ 334 oldCapacity+((oldCapacity+3)>>2)); 316 newVector>setDenseVector(vector>denseVector() + oldCapacity); 317 newVector>setIndexVector(vector>getIndices() + oldCapacity + ((oldCapacity + 3) >> 2)); 335 318 vector>checkClean(); 336 319 newVector>checkClear(); 337 320 return newVector; 338 321 } 339 static void deleteFakeVector(CoinIndexedVector * 340 CoinIndexedVector *fakeVector)341 { 342 int * 322 static void deleteFakeVector(CoinIndexedVector *vector, 323 CoinIndexedVector *fakeVector) 324 { 325 int *index = vector>getIndices(); 343 326 fakeVector>checkClear(); 344 327 fakeVector>setCapacity(0); … … 349 332 } 350 333 // Synchronize stuff 351 void 352 AbcSimplexFactorization::synchronize(const ClpFactorization * otherFactorization,const AbcSimplex * model) 353 { 354 goDenseThreshold_=otherFactorization>goDenseThreshold(); 355 goSmallThreshold_=otherFactorization>goSmallThreshold(); 356 goLongThreshold_=otherFactorization>goOslThreshold(); 334 void AbcSimplexFactorization::synchronize(const ClpFactorization *otherFactorization, const AbcSimplex *model) 335 { 336 goDenseThreshold_ = otherFactorization>goDenseThreshold(); 337 goSmallThreshold_ = otherFactorization>goSmallThreshold(); 338 goLongThreshold_ = otherFactorization>goOslThreshold(); 357 339 //forceOtherFactorization(otherFactorization>typeOfFactorization()); 358 340 goDenseOrSmall(model>numberRows()); 359 maximumPivots(static_cast< int>(otherFactorization>maximumPivots()*1.2));341 maximumPivots(static_cast< int >(otherFactorization>maximumPivots() * 1.2)); 360 342 #ifdef ABC_USE_COIN_FACTORIZATION 361 343 // redo region sizes 362 int maximumRows =model>numberRows()+maximumPivots()+1;344 int maximumRows = model>numberRows() + maximumPivots() + 1; 363 345 int currentCapacity = model>usefulArray(0)>capacity(); 364 int newCapacity = currentCapacity +maximumRows+3;365 for (int i =0;i<ABC_NUMBER_USEFUL;i++) {366 CoinPartitionedVector * 346 int newCapacity = currentCapacity + maximumRows + 3; 347 for (int i = 0; i < ABC_NUMBER_USEFUL; i++) { 348 CoinPartitionedVector *vector = model>usefulArray(i); 367 349 vector>reserve(newCapacity); 368 // zero 369 CoinZeroN(vector>getIndices(), newCapacity);350 // zero 351 CoinZeroN(vector>getIndices(), newCapacity); 370 352 // but pretend old 371 353 //vector>setCapacity(currentCapacity); … … 395 377 extern int currentTakeoutU; 396 378 #endif 397 int 398 AbcSimplexFactorization::factorize ( AbcSimplex * model, 399 int solveType, bool valuesPass) 379 int AbcSimplexFactorization::factorize(AbcSimplex *model, 380 int solveType, bool valuesPass) 400 381 { 401 382 #ifdef CLP_FACTORIZATION_INSTRUMENT 402 externalTimeStart =CoinCpuTime();403 #endif 404 model_ = model;405 AbcMatrix * 383 externalTimeStart = CoinCpuTime(); 384 #endif 385 model_ = model; 386 AbcMatrix *matrix = model>abcMatrix(); 406 387 int numberRows = model>numberRows(); 407 388 if (!numberRows) … … 411 392 #ifndef ABC_USE_COIN_FACTORIZATION 412 393 const 413 #endif 414 int * COIN_RESTRICT pivotVariable = model>pivotVariable(); 394 #endif 395 int *COIN_RESTRICT pivotVariable 396 = model>pivotVariable(); 415 397 //returns 0 okay, 1 singular, 2 too many in basis */ 416 398 // allow dense 417 int solveMode = coinAbcFactorization_>solveMode() &1;418 if (model>numberIterations() >model>baseIteration())419 solveMode = 9; 399 int solveMode = coinAbcFactorization_>solveMode() & 1; 400 if (model>numberIterations() > model>baseIteration()) 401 solveMode = 9; // was +8  this allows dense 420 402 else 421 403 solveMode = 1; // try dense … … 424 406 coinAbcFactorization_>setSolveMode(solveMode); 425 407 while (status() < 98) { 426 408 427 409 int i; 428 410 int numberBasic = 0; 429 411 // Move pivot variables across if they look good 430 int * 412 int *COIN_RESTRICT pivotTemp = model>usefulArray(0)>getIndices(); 431 413 #ifndef NDEBUG 432 414 model_>checkArrays(); 433 415 #endif 434 assert 416 assert(!model>usefulArray(0)>getNumElements()); 435 417 // Seems to prefer things in order so quickest 436 418 // way is to go though like this 437 419 for (i = 0; i < numberRows; i++) { 438 if (pivotVariable[i] <numberRows)439 420 if (pivotVariable[i] < numberRows) 421 pivotTemp[numberBasic++] = pivotVariable[i]; 440 422 } 441 423 numberSlacks_ = numberBasic; … … 443 425 */ 444 426 for (i = 0; i < numberRows; i++) { 445 if (pivotVariable[i] >=numberRows)446 427 if (pivotVariable[i] >= numberRows) 428 pivotTemp[numberBasic++] = pivotVariable[i]; 447 429 } 448 430 CoinBigIndex numberElements = numberSlacks_; 449 431 450 432 // compute how much in basis 451 433 int numberColumnBasic = numberBasic  numberSlacks_; 452 434 453 435 numberElements += matrix>countBasis(pivotTemp + numberSlacks_, 454 436 numberColumnBasic); 455 437 //printf("Basis has %d slacks  size %d\n",numberSlacks_,numberElements); 456 438 // Not needed for dense … … 463 445 coinAbcFactorization_>gutsOfInitialize(2); 464 446 #endif 465 coinAbcFactorization_>getAreas (numberRows,466 467 2 * numberElements);447 coinAbcFactorization_>getAreas(numberRows, 448 numberSlacks_ + numberColumnBasic, numberElements, 449 2 * numberElements); 468 450 #if 0 469 451 if (!model>numberIterations()) … … 472 454 // Fill in counts so we can skip part of preProcess 473 455 // This is NOT needed for dense but would be needed for later versions 474 CoinFactorizationDouble * 475 int * 476 CoinBigIndex * 477 int * 478 int * 456 CoinFactorizationDouble *COIN_RESTRICT elementU; 457 int *COIN_RESTRICT indexRowU; 458 CoinBigIndex *COIN_RESTRICT startColumnU; 459 int *COIN_RESTRICT numberInRow; 460 int *COIN_RESTRICT numberInColumn; 479 461 #define slackValue 1.0 480 462 #ifndef ABC_USE_COIN_FACTORIZATION … … 485 467 numberInColumn = coinAbcFactorization_>numberInColumn(); 486 468 coinAbcFactorization_>setNumberSlacks(numberSlacks_); 487 CoinZeroN ( numberInRow, numberRows);488 CoinZeroN ( numberInColumn, numberRows);469 CoinZeroN(numberInRow, numberRows); 470 CoinZeroN(numberInColumn, numberRows); 489 471 #else 490 472 elementU = coinAbcFactorization_>elementU(); … … 493 475 numberInRow = coinAbcFactorization_>numberInRow(); 494 476 numberInColumn = coinAbcFactorization_>numberInColumn(); 495 CoinZeroN ( numberInRow, coinAbcFactorization_>numberRows() + 1);496 CoinZeroN ( numberInColumn, coinAbcFactorization_>maximumColumnsExtra() + 1);477 CoinZeroN(numberInRow, coinAbcFactorization_>numberRows() + 1); 478 CoinZeroN(numberInColumn, coinAbcFactorization_>maximumColumnsExtra() + 1); 497 479 #endif 498 480 for (i = 0; i < numberSlacks_; i++) { … … 508 490 numberColumnBasic = numberRows  numberSlacks_; 509 491 matrix>fillBasis(pivotTemp + numberSlacks_, 510 511 512 513 514 515 516 numberElements = startColumnU[numberRows 1]517 + numberInColumn[numberRows 1];518 #ifndef ABC_USE_COIN_FACTORIZATION 519 coinAbcFactorization_>preProcess ();520 coinAbcFactorization_>factor 492 numberColumnBasic, 493 indexRowU, 494 startColumnU + numberSlacks_, 495 numberInRow, 496 numberInColumn + numberSlacks_, 497 elementU); 498 numberElements = startColumnU[numberRows  1] 499 + numberInColumn[numberRows  1]; 500 #ifndef ABC_USE_COIN_FACTORIZATION 501 coinAbcFactorization_>preProcess(); 502 coinAbcFactorization_>factor(model); 521 503 #else 522 504 // recompute number basic 523 505 numberBasic = numberSlacks_ + numberColumnBasic; 524 506 if (numberBasic) 525 numberElements = startColumnU[numberBasic 1]526 + numberInColumn[numberBasic1];507 numberElements = startColumnU[numberBasic  1] 508 + numberInColumn[numberBasic  1]; 527 509 else 528 510 numberElements = 0; 529 511 coinAbcFactorization_>setNumberElementsU(numberElements); 530 512 #ifdef CLP_FACTORIZATION_NEW_TIMING 531 lastNumberPivots_ =0;532 effectiveStartNumberU_ =numberElementsnumberRows;513 lastNumberPivots_ = 0; 514 effectiveStartNumberU_ = numberElements  numberRows; 533 515 //printf("%d slacks,%d in U at beginning\n", 534 516 //numberRowBasic,numberElements); … … 536 518 //saveFactorization("dump.d"); 537 519 if (coinAbcFactorization_>biasLU() >= 3  coinAbcFactorization_>numberRows() != coinAbcFactorization_>numberColumns()) 538 coinAbcFactorization_>preProcess ( 2);520 coinAbcFactorization_>preProcess(2); 539 521 else 540 coinAbcFactorization_>preProcess ( 3); // no row copy541 coinAbcFactorization_>factor ();522 coinAbcFactorization_>preProcess(3); // no row copy 523 coinAbcFactorization_>factor(); 542 524 #endif 543 525 #if 0 … … 548 530 } 549 531 #endif 550 if (coinAbcFactorization_>status() == 1 && 551 (coinAbcFactorization_>solveMode() & 1) != 0) { 532 if (coinAbcFactorization_>status() == 1 && (coinAbcFactorization_>solveMode() & 1) != 0) { 552 533 int solveMode = coinAbcFactorization_>solveMode(); 553 solveMode 534 solveMode; // so bottom will be 0 554 535 coinAbcFactorization_>setSolveMode(solveMode); 555 536 coinAbcFactorization_>setStatus(99); 556 537 } else if (coinAbcFactorization_>status() == 99) { 557 538 // get more memory 558 coinAbcFactorization_>areaFactor( 559 } 560 if (coinAbcFactorization_>status() == 99) 539 coinAbcFactorization_>areaFactor(coinAbcFactorization_>areaFactor() * 2.0); 540 } 541 if (coinAbcFactorization_>status() == 99) 561 542 continue; 562 543 // If we get here status is 0 or 1 … … 565 546 coinAbcFactorization_>postProcess(pivotTemp, model>pivotVariable()); 566 547 #else 567 const int * 568 const int * 548 const int *permuteBack = coinAbcFactorization_>permuteBack(); 549 const int *back = coinAbcFactorization_>pivotColumnBack(); 569 550 // Redo pivot order 570 551 for (i = 0; i < numberRows; i++) { 571 572 573 574 575 552 int k = pivotTemp[i]; 553 // so rowIsBasic[k] would be permuteBack[back[i]] 554 int j = permuteBack[back[i]]; 555 //assert (pivotVariable[j] == 1); 556 pivotVariable[j] = k; 576 557 } 577 558 // Set up permutation vector 578 559 // these arrays start off as copies of permute 579 560 // (and we could use permute_ instead of pivotColumn (not back though)) 580 ClpDisjointCopyN ( coinAbcFactorization_>permute(), numberRows , coinAbcFactorization_>pivotColumn());581 ClpDisjointCopyN ( coinAbcFactorization_>permuteBack(), numberRows , coinAbcFactorization_>pivotColumnBack());561 ClpDisjointCopyN(coinAbcFactorization_>permute(), numberRows, coinAbcFactorization_>pivotColumn()); 562 ClpDisjointCopyN(coinAbcFactorization_>permuteBack(), numberRows, coinAbcFactorization_>pivotColumnBack()); 582 563 // See if worth going sparse and when 583 564 coinAbcFactorization_>checkSparse(); 584 565 #ifdef CLP_FACTORIZATION_NEW_TIMING 585 endLengthU_ = coinAbcFactorization_>numberElements()  586 coinAbcFactorization_>numberDense()*coinAbcFactorization_>numberDense() 587 coinAbcFactorization_>numberElementsL(); 566 endLengthU_ = coinAbcFactorization_>numberElements()  coinAbcFactorization_>numberDense() * coinAbcFactorization_>numberDense() 567  coinAbcFactorization_>numberElementsL(); 588 568 #endif 589 569 #endif 590 570 model_>moveToBasic(); 591 } else if (solveType == 0  solveType == 2 /*solveType==1*/) {571 } else if (solveType == 0  solveType == 2 /*solveType==1*/) { 592 572 // Change pivotTemp to be correct list 593 573 anyChanged = true; 594 574 coinAbcFactorization_>makeNonSingular(pivotTemp); 595 const double * 596 const double * 597 double * 575 const double *COIN_RESTRICT lowerArray = model>lowerRegion(); 576 const double *COIN_RESTRICT upperArray = model>upperRegion(); 577 double *COIN_RESTRICT solution = model>solutionRegion(); 598 578 //int * pivotVariable=model_>pivotVariable(); 599 579 //int * fromExternal=model_>fromExternal(); 600 int numberTotal =model_>numberTotal();580 int numberTotal = model_>numberTotal(); 601 581 //can use external status_ 602 unsigned char * 603 CoinAbcMemset0(statusArray, numberTotal);604 for (int iRow =0;iRow<numberRows;iRow++) {605 int iPivot=pivotVariable[iRow];606 607 608 statusArray[iPivot]=1;582 unsigned char *COIN_RESTRICT statusArray = model_>statusArray(); 583 CoinAbcMemset0(statusArray, numberTotal); 584 for (int iRow = 0; iRow < numberRows; iRow++) { 585 int iPivot = pivotVariable[iRow]; 586 //if (iPivot!=pivotTemp[iRow]) 587 //printf("row %d bad pivot %d good %d\n",iRow,iPivot,pivotTemp[iRow]); 588 statusArray[iPivot] = 1; 609 589 } 610 for (int iRow =0;iRow<numberRows;iRow++) {611 int iPivot=pivotTemp[iRow];612 statusArray[iPivot]=2;590 for (int iRow = 0; iRow < numberRows; iRow++) { 591 int iPivot = pivotTemp[iRow]; 592 statusArray[iPivot] = 2; 613 593 } 614 int jPivot =0;594 int jPivot = 0; 615 595 double largeValue = model>largeValue(); 616 for (int iRow =0;iRow<numberRows;iRow++) {617 int iPivot=pivotVariable[iRow];618 if (statusArray[iPivot]==1) {619 620 while (statusArray[jPivot]!=2) {621 622 623 statusArray[iPivot]=0;624 statusArray[jPivot]=0;596 for (int iRow = 0; iRow < numberRows; iRow++) { 597 int iPivot = pivotVariable[iRow]; 598 if (statusArray[iPivot] == 1) { 599 // clean 600 while (statusArray[jPivot] != 2) { 601 jPivot++; 602 } 603 statusArray[iPivot] = 0; 604 statusArray[jPivot] = 0; 625 605 #ifndef NDEBUG 626 printf("On Row %d replacing %d by %d\n",iRow,iPivot,jPivot);627 #endif 628 629 630 631 632 633 634 if (lower!=upper) {635 636 thisStatus=AbcSimplex::atLowerBound;637 638 639 thisStatus= AbcSimplex::atUpperBound;640 641 642 643 thisStatus= AbcSimplex::isFixed;644 645 646 thisStatus=AbcSimplex::isFree;647 648 649 thisStatus=AbcSimplex::superBasic;650 651 model_>setInternalStatus(iPivot,thisStatus);652 model_>setInternalStatus(jPivot,AbcSimplex::basic);653 654 model_>swap(iRow,jPivot);655 606 printf("On Row %d replacing %d by %d\n", iRow, iPivot, jPivot); 607 #endif 608 AbcSimplex::Status thisStatus; 609 if (!valuesPass) { 610 double lower = lowerArray[iPivot]; 611 double upper = upperArray[iPivot]; 612 double value = solution[iPivot]; 613 if (lower > largeValue  upper < largeValue) { 614 if (lower != upper) { 615 if (fabs(value  lower) < fabs(value  upper)) { 616 thisStatus = AbcSimplex::atLowerBound; 617 solution[iPivot] = lower; 618 } else { 619 thisStatus = AbcSimplex::atUpperBound; 620 solution[iPivot] = upper; 621 } 622 } else { 623 thisStatus = AbcSimplex::isFixed; 624 } 625 } else { 626 thisStatus = AbcSimplex::isFree; 627 } 628 } else { 629 thisStatus = AbcSimplex::superBasic; 630 } 631 model_>setInternalStatus(iPivot, thisStatus); 632 model_>setInternalStatus(jPivot, AbcSimplex::basic); 633 // swap (solution will be wrong  but that doesn't matter as basic) 634 model_>swap(iRow, jPivot); 635 } 656 636 } 657 CoinAbcMemcpy(model_>pivotVariable(), pivotTemp,numberRows);637 CoinAbcMemcpy(model_>pivotVariable(), pivotTemp, numberRows); 658 638 #ifndef NDEBUG 659 639 model_>checkConsistentPivots(); … … 665 645 } 666 646 //coinAbcFactorization_>setSolveMode(solveMode1); 667 if ( 647 if (anyChanged && model>algorithm() < 0 && solveType > 0) { 668 648 double dummyCost; 669 static_cast< AbcSimplexDual *> (model)>changeBounds(3,dummyCost);649 static_cast< AbcSimplexDual * >(model)>changeBounds(3, dummyCost); 670 650 } 671 651 return coinAbcFactorization_>status(); … … 738 718 /* Replaces one Column to basis, 739 719 partial update already in U */ 740 void 741 AbcSimplexFactorization::replaceColumnPart3 ( const AbcSimplex * model, 742 CoinIndexedVector * regionSparse, 743 CoinIndexedVector * tableauColumn, 744 int pivotRow, 745 #ifdef ABC_LONG_FACTORIZATION 746 long 747 #endif 748 double alpha ) 720 void AbcSimplexFactorization::replaceColumnPart3(const AbcSimplex *model, 721 CoinIndexedVector *regionSparse, 722 CoinIndexedVector *tableauColumn, 723 int pivotRow, 724 #ifdef ABC_LONG_FACTORIZATION 725 long 726 #endif 727 double alpha) 749 728 { 750 729 #ifndef ABC_USE_COIN_FACTORIZATION … … 754 733 coinAbcFactorization_>setUsefulInformation(tempInfo, 1); 755 734 if (tab) 756 coinAbcFactorization_>replaceColumnPart3(model, NULL,tableauColumn,757 758 735 coinAbcFactorization_>replaceColumnPart3(model, NULL, tableauColumn, 736 pivotRow, 737 tableauColumn>denseVector()[pivotRow]); 759 738 else 760 coinAbcFactorization_>replaceColumnPart3(model, regionSparse,NULL,761 762 739 coinAbcFactorization_>replaceColumnPart3(model, regionSparse, NULL, 740 pivotRow, 741 alpha); 763 742 #else 764 743 #ifdef CLP_FACTORIZATION_NEW_TIMING 765 int nOld =0;766 int nNew =0;744 int nOld = 0; 745 int nNew = 0; 767 746 int seq; 768 const CoinPackedMatrix * matrix=model>matrix();769 const int * 770 seq =model>sequenceIn();771 if (seq >=0&&seq<model>numberColumns()+model>numberRows()) {772 if (seq <model>numberRows())773 nNew =1;747 const CoinPackedMatrix *matrix = model>matrix(); 748 const int *columnLength = matrix>getVectorLengths(); 749 seq = model>sequenceIn(); 750 if (seq >= 0 && seq < model>numberColumns() + model>numberRows()) { 751 if (seq < model>numberRows()) 752 nNew = 1; 774 753 else 775 nNew =columnLength[seqmodel>numberRows()];776 } 777 seq =model>sequenceOut();778 if (seq >=0&&seq<model>numberColumns()+model>numberRows()) {779 if (seq <model>numberRows())780 nOld =1;754 nNew = columnLength[seq  model>numberRows()]; 755 } 756 seq = model>sequenceOut(); 757 if (seq >= 0 && seq < model>numberColumns() + model>numberRows()) { 758 if (seq < model>numberRows()) 759 nOld = 1; 781 760 else 782 nOld =columnLength[seqmodel>numberRows()];783 } 784 effectiveStartNumberU_ += nNew nOld;761 nOld = columnLength[seq  model>numberRows()]; 762 } 763 effectiveStartNumberU_ += nNew  nOld; 785 764 #endif 786 765 coinAbcFactorization_>replaceColumnPart3(regionSparse, 787 788 766 pivotRow, 767 alpha); 789 768 #endif 790 769 } 791 770 /* Replaces one Column to basis, 792 771 partial update in vector */ 793 void 794 AbcSimplexFactorization::replaceColumnPart3 ( const AbcSimplex * model, 795 CoinIndexedVector * regionSparse, 796 CoinIndexedVector * tableauColumn, 797 CoinIndexedVector * partialUpdate, 798 int pivotRow, 799 #ifdef ABC_LONG_FACTORIZATION 800 long 801 #endif 802 double alpha ) 772 void AbcSimplexFactorization::replaceColumnPart3(const AbcSimplex *model, 773 CoinIndexedVector *regionSparse, 774 CoinIndexedVector *tableauColumn, 775 CoinIndexedVector *partialUpdate, 776 int pivotRow, 777 #ifdef ABC_LONG_FACTORIZATION 778 long 779 #endif 780 double alpha) 803 781 { 804 782 #ifndef ABC_USE_COIN_FACTORIZATION … … 808 786 coinAbcFactorization_>setUsefulInformation(tempInfo, 1); 809 787 if (tab) 810 coinAbcFactorization_>replaceColumnPart3(model, NULL,tableauColumn,811 812 788 coinAbcFactorization_>replaceColumnPart3(model, NULL, tableauColumn, 789 pivotRow, 790 tableauColumn>denseVector()[pivotRow]); 813 791 else 814 coinAbcFactorization_>replaceColumnPart3(model, regionSparse,NULL,partialUpdate,815 816 817 #else 818 coinAbcFactorization_>replaceColumnPart3(regionSparse,partialUpdate,819 820 792 coinAbcFactorization_>replaceColumnPart3(model, regionSparse, NULL, partialUpdate, 793 pivotRow, 794 alpha); 795 #else 796 coinAbcFactorization_>replaceColumnPart3(regionSparse, partialUpdate, 797 pivotRow, 798 alpha); 821 799 #endif 822 800 } … … 905 883 #endif 906 884 /* makes a row copy of L for speed and to allow very sparse problems */ 907 void 908 AbcSimplexFactorization::goSparse() 885 void AbcSimplexFactorization::goSparse() 909 886 { 910 887 #ifdef ABC_USE_COIN_FACTORIZATION … … 918 895 } 919 896 // Set tolerances to safer of existing and given 920 void 921 AbcSimplexFactorization::saferTolerances ( double zeroValue, 922 double pivotValue) 897 void AbcSimplexFactorization::saferTolerances(double zeroValue, 898 double pivotValue) 923 899 { 924 900 double newValue1; … … 927 903 newValue1 = zeroValue; 928 904 else 929 newValue1 = zeroTolerance() * zeroValue; 930 newValue1 = CoinMin(zeroTolerance(), newValue1);931 if (newValue1 >1.0e15)905 newValue1 = zeroTolerance() * zeroValue; 906 newValue1 = CoinMin(zeroTolerance(), newValue1); 907 if (newValue1 > 1.0e15) 932 908 zeroTolerance(newValue1); 933 double 909 double newValue2; 934 910 // better to have large tolerance even if slower 935 911 if (pivotValue > 0.0) … … 937 913 else 938 914 newValue2 = pivotTolerance() * pivotValue; 939 newValue2 = CoinMin(CoinMax(pivotTolerance(), newValue2), 0.999);940 if (newValue2 >pivotTolerance()) {915 newValue2 = CoinMin(CoinMax(pivotTolerance(), newValue2), 0.999); 916 if (newValue2 > pivotTolerance()) { 941 917 pivotTolerance(newValue2); 942 918 char line[100]; 943 sprintf(line, "new zero tolerance %g new pivot tolerance %g",944 zeroTolerance(),pivotTolerance());945 model_>messageHandler()>message(CLP_GENERAL2, *model_>messagesPointer())919 sprintf(line, "new zero tolerance %g new pivot tolerance %g", 920 zeroTolerance(), pivotTolerance()); 921 model_>messageHandler()>message(CLP_GENERAL2, *model_>messagesPointer()) 946 922 << line << CoinMessageEol; 947 923 } 948 924 } 949 925 // Sets factorization 950 void 951 AbcSimplexFactorization::setFactorization(AbcSimplexFactorization & rhs) 926 void AbcSimplexFactorization::setFactorization(AbcSimplexFactorization &rhs) 952 927 { 953 928 AbcSimplexFactorization::operator=(rhs); 954 929 } 930 931 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 932 */
Note: See TracChangeset
for help on using the changeset viewer.