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

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/AbcSimplexPrimal.cpp
r2166 r2385 108 108 Variable will change theta if currentValue  currentTheta*alpha < 0.0 109 109 */ 110 bool userChoiceValid1(const AbcSimplex * 111 112 113 114 115 110 bool userChoiceValid1(const AbcSimplex *model, 111 int sequenceOut, 112 double currentValue, 113 double currentTheta, 114 double alpha, 115 double realAlpha); 116 116 /* This returns true if chosen in/out pair valid. 117 117 The main thing to check would be variable flipping bounds may be … … 119 119 If you return false sequenceIn_ will be flagged as ineligible. 120 120 */ 121 bool userChoiceValid2(const AbcSimplex * 121 bool userChoiceValid2(const AbcSimplex *model); 122 122 /* If a good pivot then you may wish to unflag some variables. 123 123 */ 124 void userChoiceWasGood(AbcSimplex * 124 void userChoiceWasGood(AbcSimplex *model); 125 125 #endif 126 126 #ifdef TRY_SPLIT_VALUES_PASS 127 static double valuesChunk =10.0;128 static double valuesRatio =0.1;129 static int valuesStop =1;130 static int keepValuesPass =1;131 static int doOrdinaryVariables =1;127 static double valuesChunk = 10.0; 128 static double valuesRatio = 0.1; 129 static int valuesStop = 1; 130 static int keepValuesPass = 1; 131 static int doOrdinaryVariables = 1; 132 132 #endif 133 133 // primal 134 int AbcSimplexPrimal::primal (int ifValuesPass, int /*startFinishOptions*/)134 int AbcSimplexPrimal::primal(int ifValuesPass, int /*startFinishOptions*/) 135 135 { 136 136 137 137 /* 138 138 Method … … 214 214 215 215 */ 216 216 217 217 algorithm_ = +1; 218 218 moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy 219 #if ABC_PARALLEL >0219 #if ABC_PARALLEL > 0 220 220 if (parallelMode()) { 221 221 // extra regions 222 222 // for moment allow for ordered factorization 223 int length =2*numberRows_+abcFactorization_>maximumPivots();223 int length = 2 * numberRows_ + abcFactorization_>maximumPivots(); 224 224 for (int i = 0; i < 6; i++) { 225 225 delete rowArray_[i]; … … 236 236 dualTolerance_ = dblParam_[ClpDualTolerance]; 237 237 primalTolerance_ = dblParam_[ClpPrimalTolerance]; 238 currentDualTolerance_ =dualTolerance_;238 currentDualTolerance_ = dualTolerance_; 239 239 //nextCleanNonBasicIteration_=baseIteration_+numberRows_; 240 240 // Save so can see if doing after dual … … 243 243 int initialNegDjs = 1; 244 244 // copy bounds to perturbation 245 CoinAbcMemcpy(abcPerturbation_, abcLower_,numberTotal_);246 CoinAbcMemcpy(abcPerturbation_ +numberTotal_,abcUpper_,numberTotal_);245 CoinAbcMemcpy(abcPerturbation_, abcLower_, numberTotal_); 246 CoinAbcMemcpy(abcPerturbation_ + numberTotal_, abcUpper_, numberTotal_); 247 247 #if 0 248 248 if (numberRows_>80000&&numberRows_<90000) { … … 266 266 statusOfProblemInPrimal(0); 267 267 /*if (!startup(ifValuesPass))*/ { 268 268 269 269 // Set average theta 270 270 abcNonLinearCost_>setAverageTheta(1.0e3); 271 271 272 272 // Say no pivot has occurred (for steepest edge and updates) 273 273 pivotRow_ = 2; 274 274 275 275 // This says whether to restore things etc 276 276 int factorType = 0; 277 277 if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) { 278 278 perturb(0); 279 if (perturbation_ !=100) {280 281 assert(!ifValuesPass);282 283 284 285 286 287 288 289 290 291 292 293 294 279 if (perturbation_ != 100) { 280 // Can't get here if values pass 281 assert(!ifValuesPass); 282 gutsOfPrimalSolution(3); 283 if (handler_>logLevel() > 2) { 284 handler_>message(CLP_SIMPLEX_STATUS, messages_) 285 << numberIterations_ << objectiveValue(); 286 handler_>printing(sumPrimalInfeasibilities_ > 0.0) 287 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_; 288 handler_>printing(sumDualInfeasibilities_ > 0.0) 289 << sumDualInfeasibilities_ << numberDualInfeasibilities_; 290 handler_>printing(numberDualInfeasibilitiesWithoutFree_ 291 < numberDualInfeasibilities_) 292 << numberDualInfeasibilitiesWithoutFree_; 293 handler_>message() << CoinMessageEol; 294 } 295 295 } 296 296 } … … 298 298 abcProgress_.fillFromModel(this); 299 299 abcProgress_.startCheck(); 300 bool needNewWeights =false;301 double pivotTolerance =abcFactorization_>pivotTolerance();300 bool needNewWeights = false; 301 double pivotTolerance = abcFactorization_>pivotTolerance(); 302 302 /* 303 303 Status of problem: … … 316 316 // If getting nowhere  why not give it a kick 317 317 #if 1 318 if (perturbation_ < 101 && numberIterations_ == 6078/*> 2 * (numberRows_ + numberColumns_)*/ && (specialOptions_ & 4) == 0319 320 318 if (perturbation_ < 101 && numberIterations_ == 6078 /*> 2 * (numberRows_ + numberColumns_)*/ && (specialOptions_ & 4) == 0 319 && initialStatus != 10) { 320 perturb(1); 321 321 } 322 322 #endif 323 323 // If we have done no iterations  special 324 324 if (lastGoodIteration_ == numberIterations_ && factorType) 325 326 if (pivotTolerance<abcFactorization_>pivotTolerance()) {327 328 329 factorType=5;330 } 331 325 factorType = 3; 326 if (pivotTolerance < abcFactorization_>pivotTolerance()) { 327 // force factorization 328 pivotTolerance = abcFactorization_>pivotTolerance(); 329 factorType = 5; 330 } 331 332 332 // may factorize, checks if problem finished 333 333 if (factorType) 334 335 needNewWeights =false;336 if (problemStatus_ ==10 && (moreSpecialOptions_ & 2048) != 0) {337 problemStatus_=0;334 statusOfProblemInPrimal(factorType + (needNewWeights ? 10 : 0)); 335 needNewWeights = false; 336 if (problemStatus_ == 10 && (moreSpecialOptions_ & 2048) != 0) { 337 problemStatus_ = 0; 338 338 } 339 339 if (initialStatus == 10) { 340 initialStatus=1;341 342 if(initialIterations != numberIterations_) {343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 } 359 340 initialStatus = 1; 341 // cleanup phase 342 if (initialIterations != numberIterations_) { 343 if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) { 344 // getting worse  try perturbing 345 if (perturbation_ < 101 && (specialOptions_ & 4) == 0) { 346 perturb(1); 347 statusOfProblemInPrimal(factorType); 348 } 349 } 350 } else { 351 // save number of negative djs 352 if (!numberPrimalInfeasibilities_) 353 initialNegDjs = numberDualInfeasibilities_; 354 // make sure weight won't be changed 355 if (infeasibilityCost_ == 1.0e10) 356 infeasibilityCost_ = 1.000001e10; 357 } 358 } 359 360 360 // Say good factorization 361 361 factorType = 1; 362 362 363 363 // Say no pivot has occurred (for steepest edge and updates) 364 364 pivotRow_ = 2; 365 365 366 366 // exit if victory declared 367 367 if (problemStatus_ >= 0) { 368 368 #ifdef CLP_USER_DRIVEN 369 int status = 370 eventHandler_>event(ClpEventHandler::endInPrimal); 371 if (status>=0&&status<10) { 372 // carry on 373 problemStatus_=1; 374 if (status==0) 375 continue; // refactorize 376 } else if (status>=10) { 377 problemStatus_=status10; 378 break; 379 } else { 380 break; 381 } 369 int status = eventHandler_>event(ClpEventHandler::endInPrimal); 370 if (status >= 0 && status < 10) { 371 // carry on 372 problemStatus_ = 1; 373 if (status == 0) 374 continue; // refactorize 375 } else if (status >= 10) { 376 problemStatus_ = status  10; 377 break; 378 } else { 379 break; 380 } 382 381 #else 383 384 #endif 385 } 386 382 break; 383 #endif 384 } 385 387 386 // test for maximum iterations 388 387 if (hitMaximumIterations()  (ifValuesPass == 2 && firstFree_ < 0)) { 389 390 391 } 392 388 problemStatus_ = 3; 389 break; 390 } 391 393 392 if (firstFree_ < 0) { 394 395 396 397 needNewWeights=true;398 399 400 401 402 403 404 393 if (ifValuesPass) { 394 // end of values pass 395 ifValuesPass = 0; 396 needNewWeights = true; 397 int status = eventHandler_>event(ClpEventHandler::endOfValuesPass); 398 if (status >= 0) { 399 problemStatus_ = 5; 400 secondaryStatus_ = ClpEventHandler::endOfValuesPass; 401 break; 402 } 403 //#define FEB_TRY 405 404 #if 1 //def FEB_TRY 406 405 if (perturbation_ < 100 407 406 #ifdef TRY_SPLIT_VALUES_PASS 408 &&valuesStop<0409 #endif 410 411 412 #endif 413 414 } 415 if (abcNonLinearCost_>numberInfeasibilities() > 0 && !abcProgress_.initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10 407 && valuesStop < 0 408 #endif 409 ) 410 perturb(0); 411 #endif 412 } 413 } 414 if (abcNonLinearCost_>numberInfeasibilities() > 0 && !abcProgress_.initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10 416 415 #ifdef TRY_SPLIT_VALUES_PASS 417 && valuesStop<0418 #endif 419 420 421 422 CoinAbcMemcpy(djSaved_,abcDj_,numberTotal_);423 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);424 CoinAbcGatherFrom(abcCost_,costBasic_,abcPivotVariable_,numberRows_);425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 } else if(value < dualTolerance_) {448 449 450 451 452 453 454 455 456 457 } else if(value < dualTolerance_) {458 459 460 461 462 463 464 465 466 467 468 } else if(value < dualTolerance_) {469 470 471 472 473 474 475 }476 477 }416 && valuesStop < 0 417 #endif 418 ) { 419 // first time infeasible  start up weight computation 420 // compute with original costs 421 CoinAbcMemcpy(djSaved_, abcDj_, numberTotal_); 422 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 423 CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_); 424 gutsOfPrimalSolution(1); 425 int numberSame = 0; 426 int numberDifferent = 0; 427 int numberZero = 0; 428 int numberFreeSame = 0; 429 int numberFreeDifferent = 0; 430 int numberFreeZero = 0; 431 int n = 0; 432 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { 433 if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) { 434 // not basic 435 double distanceUp = abcUpper_[iSequence]  abcSolution_[iSequence]; 436 double distanceDown = abcSolution_[iSequence]  abcLower_[iSequence]; 437 double feasibleDj = abcDj_[iSequence]; 438 double infeasibleDj = djSaved_[iSequence]  feasibleDj; 439 double value = feasibleDj * infeasibleDj; 440 if (distanceUp > primalTolerance_) { 441 // Check if "free" 442 if (distanceDown > primalTolerance_) { 443 // free 444 if (value > dualTolerance_) { 445 numberFreeSame++; 446 } else if (value < dualTolerance_) { 447 numberFreeDifferent++; 448 abcDj_[n++] = feasibleDj / infeasibleDj; 449 } else { 450 numberFreeZero++; 451 } 452 } else { 453 // should not be negative 454 if (value > dualTolerance_) { 455 numberSame++; 456 } else if (value < dualTolerance_) { 457 numberDifferent++; 458 abcDj_[n++] = feasibleDj / infeasibleDj; 459 } else { 460 numberZero++; 461 } 462 } 463 } else if (distanceDown > primalTolerance_) { 464 // should not be positive 465 if (value > dualTolerance_) { 466 numberSame++; 467 } else if (value < dualTolerance_) { 468 numberDifferent++; 469 abcDj_[n++] = feasibleDj / infeasibleDj; 470 } else { 471 numberZero++; 472 } 473 } 474 } 475 abcProgress_.initialWeight_ = 1.0; 476 } 478 477 #ifdef CLP_USEFUL_PRINTOUT 479 480 numberSame,numberDifferent,numberZero,481 numberFreeSame,numberFreeDifferent,numberFreeZero);482 #endif 483 484 485 478 printf("XXXX %d same, %d different, %d zero,  free %d %d %d\n", 479 numberSame, numberDifferent, numberZero, 480 numberFreeSame, numberFreeDifferent, numberFreeZero); 481 #endif 482 // we want most to be same 483 if (n) { 484 std::sort(abcDj_, abcDj_ + n); 486 485 #ifdef CLP_USEFUL_PRINTOUT 487 488 int which = static_cast<int> ((1.0  most) * static_cast<double>(n));489 490 printf("XXXX inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take2,abcDj_[0]*infeasibilityCost_,abcDj_[n1]*infeasibilityCost_);491 #endif 492 493 494 486 double most = 0.95; 487 int which = static_cast< int >((1.0  most) * static_cast< double >(n)); 488 double take2 = abcDj_[which] * infeasibilityCost_; 489 printf("XXXX inf cost %g take %g (range %g %g)\n", infeasibilityCost_, take2, abcDj_[0] * infeasibilityCost_, abcDj_[n  1] * infeasibilityCost_); 490 #endif 491 double take = abcDj_[0] * infeasibilityCost_; 492 // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10); 493 infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10); 495 494 #ifdef CLP_USEFUL_PRINTOUT 496 printf("XXXX changing weight to %g\n",infeasibilityCost_);497 #endif 498 499 needNewWeights=true;500 501 502 495 printf("XXXX changing weight to %g\n", infeasibilityCost_); 496 #endif 497 // may need to force redo of weights 498 needNewWeights = true; 499 } 500 abcNonLinearCost_>checkInfeasibilities(0.0); 501 gutsOfPrimalSolution(3); 503 502 } 504 503 // Check event 505 504 { 506 507 508 509 510 511 505 int status = eventHandler_>event(ClpEventHandler::endOfFactorization); 506 if (status >= 0) { 507 problemStatus_ = 5; 508 secondaryStatus_ = ClpEventHandler::endOfFactorization; 509 break; 510 } 512 511 } 513 512 // Iterate … … 532 531 unflag(); 533 532 abcProgress_.clearTimesFlagged(); 534 #if ABC_PARALLEL >0533 #if ABC_PARALLEL > 0 535 534 if (parallelMode()) { 536 535 for (int i = 0; i < 6; i++) { … … 544 543 //finish(startFinishOptions); 545 544 restoreData(data); 546 setObjectiveValue(abcNonLinearCost_>feasibleReportCost() +objectiveOffset_);545 setObjectiveValue(abcNonLinearCost_>feasibleReportCost() + objectiveOffset_); 547 546 // clean up 548 if (problemStatus_ ==10)547 if (problemStatus_ == 10) 549 548 abcNonLinearCost_>checkInfeasibilities(0.0); 550 549 delete abcNonLinearCost_; 551 abcNonLinearCost_ =NULL;550 abcNonLinearCost_ = NULL; 552 551 #if 0 553 552 if (numberRows_>80000&&numberRows_<90000) { … … 575 574 +3 max iterations 576 575 */ 577 int 578 AbcSimplexPrimal::whileIterating(int valuesOption) 576 int AbcSimplexPrimal::whileIterating(int valuesOption) 579 577 { 580 578 #if 1 581 579 #define DELAYED_UPDATE 582 arrayForBtran_ =0;580 arrayForBtran_ = 0; 583 581 //setUsedArray(arrayForBtran_); 584 arrayForFtran_ =1;582 arrayForFtran_ = 1; 585 583 setUsedArray(arrayForFtran_); 586 arrayForFlipBounds_ =2;584 arrayForFlipBounds_ = 2; 587 585 setUsedArray(arrayForFlipBounds_); 588 arrayForTableauRow_ =3;586 arrayForTableauRow_ = 3; 589 587 setUsedArray(arrayForTableauRow_); 590 588 //arrayForDualColumn_=4; 591 589 //setUsedArray(arrayForDualColumn_); 592 #if ABC_PARALLEL <2593 arrayForReplaceColumn_ =4; //4;590 #if ABC_PARALLEL < 2 591 arrayForReplaceColumn_ = 4; //4; 594 592 #else 595 arrayForReplaceColumn_ =6; //4;593 arrayForReplaceColumn_ = 6; //4; 596 594 setUsedArray(arrayForReplaceColumn_); 597 595 #endif … … 605 603 if (valuesOption > 1) 606 604 superBasicType = 3; 607 int numberStartingInfeasibilities =abcNonLinearCost_>numberInfeasibilities();605 int numberStartingInfeasibilities = abcNonLinearCost_>numberInfeasibilities(); 608 606 // status stays at 1 while iterating, >=0 finished, 2 to invert 609 607 // status 3 to go to top without an invert 610 int numberFlaggedStart = abcProgress_.timesFlagged();608 int numberFlaggedStart = abcProgress_.timesFlagged(); 611 609 while (problemStatus_ == 1) { 612 610 if (!ifValuesPass) { … … 616 614 // NOTE rowArray_[0] is used by computeDuals which is a 617 615 // slow way of getting duals but might be used 618 int saveSequence =sequenceIn_;616 int saveSequence = sequenceIn_; 619 617 primalColumn(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForTableauRow_], 620 621 if (saveSequence >=0&&saveSequence!=sequenceOut_) {622 if (getInternalStatus(saveSequence)==basic)623 abcDj_[saveSequence]=0.0;624 } 625 #if ABC_NORMAL_DEBUG >0626 if (handler_>logLevel() ==63) {627 for (int i=0;i<numberTotal_;i++) {628 if (fabs(abcDj_[i])>1.0e2)629 assert (getInternalStatus(i)!=basic);630 631 double largestCost=0.0;632 double largestDj=0.0;633 double largestGoodDj=0.0;634 int iLargest=1;635 int numberInfeasibilities=0;636 double sum=0.0;637 for (int i=0;i<numberTotal_;i++) {638 if (getInternalStatus(i)==isFixed)639 640 largestCost=CoinMax(largestCost,fabs(abcCost_[i]));641 double dj=abcDj_[i];642 if (getInternalStatus(i)==atLowerBound)643 dj=dj;644 largestDj=CoinMax(largestDj,fabs(dj));645 if(largestGoodDj<dj) {646 largestGoodDj=dj;647 iLargest=i;648 649 if (getInternalStatus(i)==basic) {650 if (abcSolution_[i]>abcUpper_[i]+primalTolerance_) {651 652 sum += abcSolution_[i]abcUpper_[i]primalTolerance_;653 } else if (abcSolution_[i]<abcLower_[i]primalTolerance_) {654 655 sum = abcSolution_[i]abcLower_[i]+primalTolerance_;656 657 658 659 660 661 if (iLargest>=numberRows_) {662 xx='C';663 kLargest=iLargestnumberRows_;664 665 xx='R';666 kLargest=iLargest;667 668 669 largestCost,largestDj,largestGoodDj,iLargest,xx,kLargest,numberInfeasibilities,sum,670 671 assert (getInternalStatus(iLargest)!=basic);618 &usefulArray_[arrayForFlipBounds_]); 619 if (saveSequence >= 0 && saveSequence != sequenceOut_) { 620 if (getInternalStatus(saveSequence) == basic) 621 abcDj_[saveSequence] = 0.0; 622 } 623 #if ABC_NORMAL_DEBUG > 0 624 if (handler_>logLevel() == 63) { 625 for (int i = 0; i < numberTotal_; i++) { 626 if (fabs(abcDj_[i]) > 1.0e2) 627 assert(getInternalStatus(i) != basic); 628 } 629 double largestCost = 0.0; 630 double largestDj = 0.0; 631 double largestGoodDj = 0.0; 632 int iLargest = 1; 633 int numberInfeasibilities = 0; 634 double sum = 0.0; 635 for (int i = 0; i < numberTotal_; i++) { 636 if (getInternalStatus(i) == isFixed) 637 continue; 638 largestCost = CoinMax(largestCost, fabs(abcCost_[i])); 639 double dj = abcDj_[i]; 640 if (getInternalStatus(i) == atLowerBound) 641 dj = dj; 642 largestDj = CoinMax(largestDj, fabs(dj)); 643 if (largestGoodDj < dj) { 644 largestGoodDj = dj; 645 iLargest = i; 646 } 647 if (getInternalStatus(i) == basic) { 648 if (abcSolution_[i] > abcUpper_[i] + primalTolerance_) { 649 numberInfeasibilities++; 650 sum += abcSolution_[i]  abcUpper_[i]  primalTolerance_; 651 } else if (abcSolution_[i] < abcLower_[i]  primalTolerance_) { 652 numberInfeasibilities++; 653 sum = abcSolution_[i]  abcLower_[i] + primalTolerance_; 654 } 655 } 656 } 657 char xx; 658 int kLargest; 659 if (iLargest >= numberRows_) { 660 xx = 'C'; 661 kLargest = iLargest  numberRows_; 662 } else { 663 xx = 'R'; 664 kLargest = iLargest; 665 } 666 printf("largest cost %g, largest dj %g best %g (%d==%c%d)  %d infeasible (sum %g) nonlininf %d\n", 667 largestCost, largestDj, largestGoodDj, iLargest, xx, kLargest, numberInfeasibilities, sum, 668 abcNonLinearCost_>numberInfeasibilities()); 669 assert(getInternalStatus(iLargest) != basic); 672 670 } 673 671 #endif 674 672 } else { 675 673 // in values pass 676 for (int i =0;i<4;i++)677 multipleSequenceIn_[i]=1;674 for (int i = 0; i < 4; i++) 675 multipleSequenceIn_[i] = 1; 678 676 int sequenceIn = nextSuperBasic(superBasicType, &usefulArray_[arrayForFlipBounds_]); 679 677 if (valuesOption > 1) 680 678 superBasicType = 2; 681 679 if (sequenceIn < 0) { 682 683 sequenceIn_=1;684 685 686 680 // end of values pass  initialize weights etc 681 sequenceIn_ = 1; 682 handler_>message(CLP_END_VALUES_PASS, messages_) 683 << numberIterations_; 684 stateOfProblem_ &= ~VALUES_PASS; 687 685 #ifdef TRY_SPLIT_VALUES_PASS 688 valuesStop=numberIterations_+doOrdinaryVariables;689 #endif 690 691 692 693 694 695 696 697 698 699 700 701 686 valuesStop = numberIterations_ + doOrdinaryVariables; 687 #endif 688 abcPrimalColumnPivot_>saveWeights(this, 5); 689 problemStatus_ = 2; // factorize now 690 pivotRow_ = 1; // say no weights update 691 returnCode = 4; 692 // Clean up 693 for (int i = 0; i < numberTotal_; i++) { 694 if (getInternalStatus(i) == atLowerBound  getInternalStatus(i) == isFixed) 695 abcSolution_[i] = abcLower_[i]; 696 else if (getInternalStatus(i) == atUpperBound) 697 abcSolution_[i] = abcUpper_[i]; 698 } 699 break; 702 700 } else { 703 704 705 706 707 708 709 710 if (maximumIterations()==100000) {711 multipleSequenceIn_[0]=sequenceIn;712 for (int i=1;i<4;i++) {713 714 if (sequenceIn>=0)715 multipleSequenceIn_[i]=sequenceIn;716 717 718 719 701 // normal 702 sequenceIn_ = sequenceIn; 703 valueIn_ = abcSolution_[sequenceIn_]; 704 lowerIn_ = abcLower_[sequenceIn_]; 705 upperIn_ = abcUpper_[sequenceIn_]; 706 dualIn_ = abcDj_[sequenceIn_]; 707 // see if any more 708 if (maximumIterations() == 100000) { 709 multipleSequenceIn_[0] = sequenceIn; 710 for (int i = 1; i < 4; i++) { 711 int sequenceIn = nextSuperBasic(superBasicType, &usefulArray_[arrayForFlipBounds_]); 712 if (sequenceIn >= 0) 713 multipleSequenceIn_[i] = sequenceIn; 714 else 715 break; 716 } 717 } 720 718 } 721 719 } … … 725 723 if (sequenceIn_ >= 0) { 726 724 // we found a pivot column 727 assert 725 assert(!flagged(sequenceIn_)); 728 726 //#define MULTIPLE_PRICE 729 727 // do second half of iteration 730 if (multipleSequenceIn_[1] ==1maximumIterations()!=100000) {731 728 if (multipleSequenceIn_[1] == 1  maximumIterations() != 100000) { 729 returnCode = pivotResult(ifValuesPass); 732 730 } else { 733 if (multipleSequenceIn_[0]<0)734 multipleSequenceIn_[0]=sequenceIn_;735 731 if (multipleSequenceIn_[0] < 0) 732 multipleSequenceIn_[0] = sequenceIn_; 733 returnCode = pivotResult4(ifValuesPass); 736 734 #ifdef MULTIPLE_PRICE 737 if (sequenceIn_>=0)738 739 #endif 740 } 741 if (numberStartingInfeasibilities&&!abcNonLinearCost_>numberInfeasibilities()) {742 743 if (abcFactorization_>pivots()>2&&numberIterations_>(numberRows_+numberColumns_)&&(stateOfProblem_&PESSIMISTIC)!=0)744 returnCode=2; // refactorize  maybe just after n iterations735 if (sequenceIn_ >= 0) 736 returnCode = pivotResult(ifValuesPass); 737 #endif 738 } 739 if (numberStartingInfeasibilities && !abcNonLinearCost_>numberInfeasibilities()) { 740 //if (abcFactorization_>pivots()>200&&numberIterations_>2*(numberRows_+numberColumns_)) 741 if (abcFactorization_>pivots() > 2 && numberIterations_ > (numberRows_ + numberColumns_) && (stateOfProblem_ & PESSIMISTIC) != 0) 742 returnCode = 2; // refactorize  maybe just after n iterations 745 743 } 746 744 if (returnCode < 1 && returnCode > 5) { 747 745 problemStatus_ = 2; // 748 746 } else if (returnCode == 5) { 749 if (abcProgress_.timesFlagged()>10+numberFlaggedStart)750 problemStatus_ =2;751 752 753 754 755 747 if (abcProgress_.timesFlagged() > 10 + numberFlaggedStart) 748 problemStatus_ = 2; 749 if ((moreSpecialOptions_ & 16) == 0 && abcFactorization_>pivots()) { 750 moreSpecialOptions_ = 16; 751 problemStatus_ = 2; 752 } 753 // otherwise something flagged  continue; 756 754 } else if (returnCode == 2) { 757 755 problemStatus_ = 5; // looks unbounded 758 756 } else if (returnCode == 4) { 759 757 problemStatus_ = 2; // looks unbounded but has iterated 760 758 } else if (returnCode != 1) { 761 762 763 759 assert(returnCode == 3); 760 if (problemStatus_ != 5) 761 problemStatus_ = 3; 764 762 } 765 763 } else { 766 764 // no pivot column 767 #if ABC_NORMAL_DEBUG >3765 #if ABC_NORMAL_DEBUG > 3 768 766 if (handler_>logLevel() & 32) 769 767 printf("** no column pivot\n"); 770 768 #endif 771 769 if (abcNonLinearCost_>numberInfeasibilities()) 772 770 problemStatus_ = 4; // might be infeasible 773 771 // Force to refactorize early next time 774 772 int numberPivots = abcFactorization_>pivots(); … … 776 774 #ifdef CLP_USER_DRIVEN 777 775 // If large number of pivots trap later? 778 if (problemStatus_ ==1 && numberPivots<13000) {779 780 if (status>=0&&status<10) {781 782 problemStatus_=1;783 if (status==0) 784 785 } else if (status>=10) {786 problemStatus_=status10;787 788 789 790 791 776 if (problemStatus_ == 1 && numberPivots < 13000) { 777 int status = eventHandler_>event(ClpEventHandler::noCandidateInPrimal); 778 if (status >= 0 && status < 10) { 779 // carry on 780 problemStatus_ = 1; 781 if (status == 0) 782 break; 783 } else if (status >= 10) { 784 problemStatus_ = status  10; 785 break; 786 } else { 787 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); 788 break; 789 } 792 790 } 793 791 #else … … 802 800 } 803 801 /* Checks if finished. Updates status */ 804 void 805 AbcSimplexPrimal::statusOfProblemInPrimal(int type) 802 void AbcSimplexPrimal::statusOfProblemInPrimal(int type) 806 803 { 807 804 int saveFirstFree = firstFree_; … … 813 810 #endif 814 811 //int saveType=type; 815 if (type >3) {812 if (type > 3) { 816 813 // force factorization 817 numberPivots =9999999;818 if (type >10)819 type = 10;820 else 814 numberPivots = 9999999; 815 if (type > 10) 816 type = 10; 817 else 821 818 type &= 3; 822 819 } 823 820 #ifndef TRY_ABC_GUS 824 bool doFactorization = (type != 3&&(numberPivots>dontFactorizePivots_numberIterations_==baseIteration_problemStatus_==10));821 bool doFactorization = (type != 3 && (numberPivots > dontFactorizePivots_  numberIterations_ == baseIteration_  problemStatus_ == 10)); 825 822 #else 826 bool doFactorization = (type != 3&&(numberPivots>dontFactorizePivots_numberIterations_==baseIteration_));827 #endif 828 bool doWeights = doFactorization problemStatus_==10;823 bool doFactorization = (type != 3 && (numberPivots > dontFactorizePivots_  numberIterations_ == baseIteration_)); 824 #endif 825 bool doWeights = doFactorization  problemStatus_ == 10; 829 826 if (type == 2) { 830 827 // trouble  restore solution … … 841 838 #endif 842 839 #if CLP_CAUTION 843 double lastAverageInfeasibility = sumDualInfeasibilities_ / 844 static_cast<double>(numberDualInfeasibilities_ + 1); 845 #endif 846 if (numberIterations_&&type) { 840 double lastAverageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 1); 841 #endif 842 if (numberIterations_ && type) { 847 843 lastSumInfeasibility = abcNonLinearCost_>sumInfeasibilities(); 848 844 lastNumberInfeasibility = abcNonLinearCost_>numberInfeasibilities(); 849 845 } else { 850 lastAverageInfeasibility =1.0e10;851 } 852 bool ifValuesPass =(stateOfProblem_&VALUES_PASS)!=0;853 bool takenAction =false;854 double sumInfeasibility =0.0;846 lastAverageInfeasibility = 1.0e10; 847 } 848 bool ifValuesPass = (stateOfProblem_ & VALUES_PASS) != 0; 849 bool takenAction = false; 850 double sumInfeasibility = 0.0; 855 851 if (problemStatus_ > 3  problemStatus_ == 4) { 856 852 // factorize … … 861 857 if (doFactorization) 862 858 abcPrimalColumnPivot_>saveWeights(this, 1); 863 859 864 860 if (!type) { 865 861 // be optimistic … … 871 867 if (doFactorization) { 872 868 // is factorization okay? 873 int solveType = ifValuesPass ? 11 : 1;869 int solveType = ifValuesPass ? 11 : 1; 874 870 if (!type) 875 871 solveType++; 876 872 int factorStatus = internalFactorize(solveType); 877 873 if (factorStatus) { 878 879 880 #if ABC_NORMAL_DEBUG >0881 882 #endif 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 874 if (type != 1  largestPrimalError_ > 1.0e3 875  largestDualError_ > 1.0e3) { 876 #if ABC_NORMAL_DEBUG > 0 877 printf("Bad initial basis\n"); 878 #endif 879 internalFactorize(2); 880 } else { 881 // no  restore previous basis 882 // Keep any flagged variables 883 for (int i = 0; i < numberTotal_; i++) { 884 if (flagged(i)) 885 internalStatusSaved_[i] = 64; //say flagged 886 } 887 restoreGoodStatus(1); 888 if (numberPivots <= 1) { 889 // throw out something 890 if (sequenceIn_ >= 0 && getInternalStatus(sequenceIn_) != basic) { 891 setFlagged(sequenceIn_); 892 } else if (sequenceOut_ >= 0 && getInternalStatus(sequenceOut_) != basic) { 893 setFlagged(sequenceOut_); 894 } 895 abcProgress_.incrementTimesFlagged(); 896 double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), abcFactorization_>pivotTolerance()); 897 abcFactorization_>pivotTolerance(newTolerance); 898 } else { 899 // Go to safe 900 abcFactorization_>pivotTolerance(0.99); 901 } 902 forceFactorization_ = 1; // a bit drastic but .. 903 type = 2; 904 abcNonLinearCost_>checkInfeasibilities(); 905 if (internalFactorize(2) != 0) { 906 largestPrimalError_ = 1.0e4; // force other type 907 } 908 } 909 changeMade_++; // say change made 914 910 } 915 911 } 916 912 if (problemStatus_ != 4) 917 913 problemStatus_ = 3; 918 if(abcProgress_.realInfeasibility_[0]<1.0e1 && 919 primalTolerance_==1.0e7&&abcProgress_.iterationNumber_[0]>0&& 920 abcProgress_.iterationNumber_[CLP_PROGRESS1]abcProgress_.iterationNumber_[0]>25) { 914 if (abcProgress_.realInfeasibility_[0] < 1.0e1 && primalTolerance_ == 1.0e7 && abcProgress_.iterationNumber_[0] > 0 && abcProgress_.iterationNumber_[CLP_PROGRESS  1]  abcProgress_.iterationNumber_[0] > 25) { 921 915 // default  so user did not set 922 916 int iP; 923 double minAverage =COIN_DBL_MAX;924 double maxAverage =0.0;925 for (iP =0;iP<CLP_PROGRESS;iP++) {926 int n=abcProgress_.numberInfeasibilities_[iP];927 928 929 930 double average=abcProgress_.realInfeasibility_[iP];931 if (average>0.1)932 933 average /= static_cast<double>(n);934 minAverage=CoinMin(minAverage,average);935 maxAverage=CoinMax(maxAverage,average);936 937 } 938 if (iP ==CLP_PROGRESS&&minAverage<1.0e5&&maxAverage<1.0e3) {939 940 #if CBC_USEFUL_PRINTING >0941 942 #endif 943 primalTolerance_=1.0e6;944 dblParam_[ClpPrimalTolerance]=1.0e6;945 917 double minAverage = COIN_DBL_MAX; 918 double maxAverage = 0.0; 919 for (iP = 0; iP < CLP_PROGRESS; iP++) { 920 int n = abcProgress_.numberInfeasibilities_[iP]; 921 if (!n) { 922 break; 923 } else { 924 double average = abcProgress_.realInfeasibility_[iP]; 925 if (average > 0.1) 926 break; 927 average /= static_cast< double >(n); 928 minAverage = CoinMin(minAverage, average); 929 maxAverage = CoinMax(maxAverage, average); 930 } 931 } 932 if (iP == CLP_PROGRESS && minAverage < 1.0e5 && maxAverage < 1.0e3) { 933 // change tolerance 934 #if CBC_USEFUL_PRINTING > 0 935 printf("CCchanging tolerance\n"); 936 #endif 937 primalTolerance_ = 1.0e6; 938 dblParam_[ClpPrimalTolerance] = 1.0e6; 939 moreSpecialOptions_ = 4194304; 946 940 } 947 941 } … … 950 944 // put back original costs and then check 951 945 // createRim(4); // costs do not change 952 if (ifValuesPass &&numberIterations_==baseIteration_) {946 if (ifValuesPass && numberIterations_ == baseIteration_) { 953 947 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 954 948 lastSumInfeasibility = abcNonLinearCost_>largestInfeasibility(); … … 960 954 } 961 955 #endif 962 if (ifValuesPass &&numberIterations_==baseIteration_) {963 double * 956 if (ifValuesPass && numberIterations_ == baseIteration_) { 957 double *save = new double[numberRows_]; 964 958 int numberOut = 1; 965 959 while (numberOut) { 966 for (int iRow = 0; iRow < numberRows_; iRow++) { 967 int iPivot = abcPivotVariable_[iRow]; 968 save[iRow] = abcSolution_[iPivot]; 969 } 970 gutsOfPrimalSolution(2); 971 double badInfeasibility = abcNonLinearCost_>largestInfeasibility(); 972 numberOut = 0; 973 // But may be very large rhs etc 974 double useError = CoinMin(largestPrimalError_, 975 1.0e5 / CoinAbcMaximumAbsElement(abcSolution_, numberTotal_)); 976 if ((lastSumInfeasibility < incomingInfeasibility_  badInfeasibility > 977 (CoinMax(10.0, 100.0 * lastSumInfeasibility))) 978 && (badInfeasibility > 10.0 useError > 1.0e3)) { 979 //printf("Original largest infeas %g, now %g, primalError %g\n", 980 // lastSumInfeasibility,abcNonLinearCost_>largestInfeasibility(), 981 // largestPrimalError_); 982 // throw out up to 1000 structurals 983 int * sort = new int[numberRows_]; 984 // first put back solution and store difference 985 for (int iRow = 0; iRow < numberRows_; iRow++) { 986 int iPivot = abcPivotVariable_[iRow]; 987 double difference = fabs(abcSolution_[iPivot]  save[iRow]); 988 abcSolution_[iPivot] = save[iRow]; 989 save[iRow] = difference; 990 } 991 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 992 if (handler_>logLevel()>1) 993 printf("Largest infeasibility %g\n",abcNonLinearCost_>largestInfeasibility()); 994 int numberBasic = 0; 995 for (int iRow = 0; iRow < numberRows_; iRow++) { 996 int iPivot = abcPivotVariable_[iRow]; 997 if (iPivot >= numberRows_) { 998 // column 999 double difference = save[iRow]; 1000 if (difference > 1.0e4) { 1001 sort[numberOut] = iRow; 1002 save[numberOut++] = difference; 1003 if (getInternalStatus(iPivot) == basic) 1004 numberBasic++; 1005 } 1006 } 1007 } 1008 if (!numberBasic) { 1009 //printf("no errors on basic  going to all slack  numberOut %d\n",numberOut); 1010 // allow 1011 numberOut = 0; 1012 } 1013 CoinSort_2(save, save + numberOut, sort); 1014 numberOut = CoinMin(1000, numberOut); 1015 // for now bring in any slack 1016 int jRow=0; 1017 for (int i = 0; i < numberOut; i++) { 1018 int iRow = sort[i]; 1019 int iColumn = abcPivotVariable_[iRow]; 1020 assert (getInternalStatus(iColumn)==basic); 1021 setInternalStatus(iColumn, superBasic); 1022 while (getInternalStatus(jRow)==basic) 1023 jRow++; 1024 setInternalStatus(jRow, basic); 1025 abcPivotVariable_[iRow] = jRow; 1026 if (fabs(abcSolution_[iColumn]) > 1.0e10) { 1027 if (abcUpper_[iColumn] < 0.0) { 1028 abcSolution_[iColumn] = abcUpper_[iColumn]; 1029 } else if (abcLower_[iColumn] > 0.0) { 1030 abcSolution_[iColumn] = abcLower_[iColumn]; 1031 } else { 1032 abcSolution_[iColumn] = 0.0; 1033 } 1034 } 1035 } 1036 delete [] sort; 1037 } 1038 if (numberOut) { 1039 double savePivotTolerance = abcFactorization_>pivotTolerance(); 1040 abcFactorization_>pivotTolerance(0.99); 1041 int factorStatus = internalFactorize(12); 1042 abcFactorization_>pivotTolerance(savePivotTolerance); 1043 assert (!factorStatus); 1044 } 1045 } 1046 delete [] save; 960 for (int iRow = 0; iRow < numberRows_; iRow++) { 961 int iPivot = abcPivotVariable_[iRow]; 962 save[iRow] = abcSolution_[iPivot]; 963 } 964 gutsOfPrimalSolution(2); 965 double badInfeasibility = abcNonLinearCost_>largestInfeasibility(); 966 numberOut = 0; 967 // But may be very large rhs etc 968 double useError = CoinMin(largestPrimalError_, 969 1.0e5 / CoinAbcMaximumAbsElement(abcSolution_, numberTotal_)); 970 if ((lastSumInfeasibility < incomingInfeasibility_  badInfeasibility > (CoinMax(10.0, 100.0 * lastSumInfeasibility))) 971 && (badInfeasibility > 10.0  useError > 1.0e3)) { 972 //printf("Original largest infeas %g, now %g, primalError %g\n", 973 // lastSumInfeasibility,abcNonLinearCost_>largestInfeasibility(), 974 // largestPrimalError_); 975 // throw out up to 1000 structurals 976 int *sort = new int[numberRows_]; 977 // first put back solution and store difference 978 for (int iRow = 0; iRow < numberRows_; iRow++) { 979 int iPivot = abcPivotVariable_[iRow]; 980 double difference = fabs(abcSolution_[iPivot]  save[iRow]); 981 abcSolution_[iPivot] = save[iRow]; 982 save[iRow] = difference; 983 } 984 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 985 if (handler_>logLevel() > 1) 986 printf("Largest infeasibility %g\n", abcNonLinearCost_>largestInfeasibility()); 987 int numberBasic = 0; 988 for (int iRow = 0; iRow < numberRows_; iRow++) { 989 int iPivot = abcPivotVariable_[iRow]; 990 if (iPivot >= numberRows_) { 991 // column 992 double difference = save[iRow]; 993 if (difference > 1.0e4) { 994 sort[numberOut] = iRow; 995 save[numberOut++] = difference; 996 if (getInternalStatus(iPivot) == basic) 997 numberBasic++; 998 } 999 } 1000 } 1001 if (!numberBasic) { 1002 //printf("no errors on basic  going to all slack  numberOut %d\n",numberOut); 1003 // allow 1004 numberOut = 0; 1005 } 1006 CoinSort_2(save, save + numberOut, sort); 1007 numberOut = CoinMin(1000, numberOut); 1008 // for now bring in any slack 1009 int jRow = 0; 1010 for (int i = 0; i < numberOut; i++) { 1011 int iRow = sort[i]; 1012 int iColumn = abcPivotVariable_[iRow]; 1013 assert(getInternalStatus(iColumn) == basic); 1014 setInternalStatus(iColumn, superBasic); 1015 while (getInternalStatus(jRow) == basic) 1016 jRow++; 1017 setInternalStatus(jRow, basic); 1018 abcPivotVariable_[iRow] = jRow; 1019 if (fabs(abcSolution_[iColumn]) > 1.0e10) { 1020 if (abcUpper_[iColumn] < 0.0) { 1021 abcSolution_[iColumn] = abcUpper_[iColumn]; 1022 } else if (abcLower_[iColumn] > 0.0) { 1023 abcSolution_[iColumn] = abcLower_[iColumn]; 1024 } else { 1025 abcSolution_[iColumn] = 0.0; 1026 } 1027 } 1028 } 1029 delete[] sort; 1030 } 1031 if (numberOut) { 1032 double savePivotTolerance = abcFactorization_>pivotTolerance(); 1033 abcFactorization_>pivotTolerance(0.99); 1034 int factorStatus = internalFactorize(12); 1035 abcFactorization_>pivotTolerance(savePivotTolerance); 1036 assert(!factorStatus); 1037 } 1038 } 1039 delete[] save; 1047 1040 } 1048 1041 gutsOfPrimalSolution(3); 1049 sumInfeasibility = 1042 sumInfeasibility = abcNonLinearCost_>sumInfeasibilities(); 1050 1043 int reason2 = 0; 1051 1044 #if CLP_CAUTION 1052 #if CLP_CAUTION ==21045 #if CLP_CAUTION == 2 1053 1046 double test2 = 1.0e5; 1054 1047 #else 1055 1048 double test2 = 1.0e1; 1056 1049 #endif 1057 if (!lastSumInfeasibility && sumInfeasibility>1.0e3*primalTolerance_ && 1058 lastAverageInfeasibility < test2 && numberPivots > 10&&!ifValuesPass) 1050 if (!lastSumInfeasibility && sumInfeasibility > 1.0e3 * primalTolerance_ && lastAverageInfeasibility < test2 && numberPivots > 10 && !ifValuesPass) 1059 1051 reason2 = 3; 1060 if (lastSumInfeasibility < 1.0e6 && sumInfeasibility > 1.0e3 && 1061 numberPivots > 10&&!ifValuesPass) 1052 if (lastSumInfeasibility < 1.0e6 && sumInfeasibility > 1.0e3 && numberPivots > 10 && !ifValuesPass) 1062 1053 reason2 = 4; 1063 1054 #endif 1064 1055 if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility 1065 && abcFactorization_>pivotTolerance() < 0.11)  1066 1056 && abcFactorization_>pivotTolerance() < 0.11) 1057  (largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10)) 1067 1058 reason2 = 2; 1068 1059 if (reason2) { 1069 takenAction =true;1060 takenAction = true; 1070 1061 problemStatus_ = tentativeStatus; 1071 1062 doFactorization = true; 1072 1063 if (numberPivots) { 1073 1074 1075 1076 sequenceIn_=1;1077 sequenceOut_=1;1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 sumInfeasibility =abcNonLinearCost_>sumInfeasibilities();1064 // go back 1065 // trouble  restore solution 1066 restoreGoodStatus(1); 1067 sequenceIn_ = 1; 1068 sequenceOut_ = 1; 1069 abcNonLinearCost_>checkInfeasibilities(); 1070 if (reason2 < 3) { 1071 // Go to safe 1072 abcFactorization_>pivotTolerance(CoinMin(0.99, 1.01 * abcFactorization_>pivotTolerance())); 1073 forceFactorization_ = 1; // a bit drastic but .. 1074 } else if (forceFactorization_ < 0) { 1075 forceFactorization_ = CoinMin(numberPivots / 4, 100); 1076 } else { 1077 forceFactorization_ = CoinMin(forceFactorization_, 1078 CoinMax(3, numberPivots / 4)); 1079 } 1080 pivotRow_ = 1; // say no weights update 1081 changeMade_++; // say change made 1082 if (numberPivots == 1) { 1083 // throw out something 1084 if (sequenceIn_ >= 0 && getInternalStatus(sequenceIn_) != basic) { 1085 setFlagged(sequenceIn_); 1086 } else if (sequenceOut_ >= 0 && getInternalStatus(sequenceOut_) != basic) { 1087 setFlagged(sequenceOut_); 1088 } 1089 abcProgress_.incrementTimesFlagged(); 1090 } 1091 type = 2; // so will restore weights 1092 if (internalFactorize(2) != 0) { 1093 largestPrimalError_ = 1.0e4; // force other type 1094 } 1095 numberPivots = 0; 1096 gutsOfPrimalSolution(3); 1097 sumInfeasibility = abcNonLinearCost_>sumInfeasibilities(); 1107 1098 } 1108 1099 } … … 1117 1108 // If in primal and small dj give up 1118 1109 if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) { 1119 double average = sumDualInfeasibilities_ / (static_cast< double>(numberDualInfeasibilities_));1110 double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_)); 1120 1111 if (numberIterations_ > 300 && average < 1.0e4) { 1121 1112 numberDualInfeasibilities_ = 0; … … 1125 1116 // Check if looping 1126 1117 int loop; 1127 ifValuesPass =firstFree_>=0;1118 ifValuesPass = firstFree_ >= 0; 1128 1119 if (type != 2 && !ifValuesPass) 1129 1120 loop = abcProgress_.looping(); … … 1131 1122 loop = 1; 1132 1123 if ((moreSpecialOptions_ & 2048) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) { 1133 double average = sumDualInfeasibilities_ / (static_cast<double> (numberDualInfeasibilities_)); 1134 if (abcProgress_.lastIterationNumber(2)==numberIterations_&& 1135 ((abcProgress_.timesFlagged()>2&&average < 1.0e1) 1136 abcProgress_.timesFlagged()>20)) { 1124 double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_)); 1125 if (abcProgress_.lastIterationNumber(2) == numberIterations_ && ((abcProgress_.timesFlagged() > 2 && average < 1.0e1)  abcProgress_.timesFlagged() > 20)) { 1137 1126 numberDualInfeasibilities_ = 0; 1138 1127 sumDualInfeasibilities_ = 0.0; 1139 problemStatus_ =3;1140 loop =0;1128 problemStatus_ = 3; 1129 loop = 0; 1141 1130 } 1142 1131 } … … 1154 1143 } 1155 1144 problemStatus_ = 10; // instead  try other algorithm 1156 return 1145 return; 1157 1146 } else if (loop < 1) { 1158 1147 // Is it time for drastic measures 1159 if (abcNonLinearCost_>numberInfeasibilities() && abcProgress_.badTimes() > 5 && 1160 abcProgress_.oddState() < 10 && abcProgress_.oddState() >= 0) { 1148 if (abcNonLinearCost_>numberInfeasibilities() && abcProgress_.badTimes() > 5 && abcProgress_.oddState() < 10 && abcProgress_.oddState() >= 0) { 1161 1149 abcProgress_.newOddState(); 1162 1150 abcNonLinearCost_>zapCosts(); … … 1174 1162 } 1175 1163 abcProgress_.modifyObjective(abcNonLinearCost_>feasibleReportCost()); 1176 if (!lastNumberInfeasibility && sumInfeasibility && numberPivots > 1&&!ifValuesPass&& 1177 !takenAction&& 1178 abcProgress_.lastObjective(0)>=abcProgress_.lastObjective(CLP_PROGRESS1)) { 1164 if (!lastNumberInfeasibility && sumInfeasibility && numberPivots > 1 && !ifValuesPass && !takenAction && abcProgress_.lastObjective(0) >= abcProgress_.lastObjective(CLP_PROGRESS  1)) { 1179 1165 // first increase minimumThetaMovement_; 1180 1166 // be more careful 1181 1167 //abcFactorization_>saferTolerances (0.99, 1.002); 1182 #if ABC_NORMAL_DEBUG >01168 #if ABC_NORMAL_DEBUG > 0 1183 1169 if (handler_>logLevel() == 63) 1184 printf("thought feasible but sum is %g force %d minimum theta %g\n", sumInfeasibility, 1185 forceFactorization_,minimumThetaMovement_);1170 printf("thought feasible but sum is %g force %d minimum theta %g\n", sumInfeasibility, 1171 forceFactorization_, minimumThetaMovement_); 1186 1172 #endif 1187 1173 stateOfProblem_ = PESSIMISTIC; 1188 if (minimumThetaMovement_ <1.0e10) {1174 if (minimumThetaMovement_ < 1.0e10) { 1189 1175 minimumThetaMovement_ *= 1.1; 1190 1176 } else if (0) { 1191 1177 int maxFactor = abcFactorization_>maximumPivots(); 1192 1178 if (maxFactor > 10) { 1193 1194 1195 1196 1197 1179 if (forceFactorization_ < 0) 1180 forceFactorization_ = maxFactor; 1181 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2)); 1182 if (handler_>logLevel() == 63) 1183 printf("Reducing factorization frequency to %d\n", forceFactorization_); 1198 1184 } 1199 1185 } … … 1205 1191 //problemStatus_=1;; 1206 1192 progressFlag_ = 0; //reset progress flag 1207 1193 1208 1194 handler_>message(CLP_SIMPLEX_STATUS, messages_) 1209 << numberIterations_ << abcNonLinearCost_>feasibleReportCost() +objectiveOffset_;1195 << numberIterations_ << abcNonLinearCost_>feasibleReportCost() + objectiveOffset_; 1210 1196 handler_>printing(abcNonLinearCost_>numberInfeasibilities() > 0) 1211 1197 << abcNonLinearCost_>sumInfeasibilities() << abcNonLinearCost_>numberInfeasibilities(); … … 1213 1199 << sumDualInfeasibilities_ << numberDualInfeasibilities_; 1214 1200 handler_>printing(numberDualInfeasibilitiesWithoutFree_ 1215 1201 < numberDualInfeasibilities_) 1216 1202 << numberDualInfeasibilitiesWithoutFree_; 1217 1203 handler_>message() << CoinMessageEol; … … 1244 1230 lastObj3 += infeasibilityCost_ * 2.0 * lastInf3; 1245 1231 if (lastObj < thisObj  1.0e5 * CoinMax(fabs(thisObj), fabs(lastObj))  1.0e7 1246 1247 #if ABC_NORMAL_DEBUG >01232 && firstFree_ < 0 && thisInf >= lastInf) { 1233 #if ABC_NORMAL_DEBUG > 0 1248 1234 if (handler_>logLevel() == 63) 1249 1235 printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_); 1250 1236 #endif 1251 1237 int maxFactor = abcFactorization_>maximumPivots(); 1252 1238 if (maxFactor > 10) { 1253 1254 1255 1256 #if ABC_NORMAL_DEBUG >01257 1258 1239 if (forceFactorization_ < 0) 1240 forceFactorization_ = maxFactor; 1241 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2)); 1242 #if ABC_NORMAL_DEBUG > 0 1243 if (handler_>logLevel() == 63) 1244 printf("Reducing factorization frequency to %d\n", forceFactorization_); 1259 1245 #endif 1260 1246 } 1261 1247 } else if (lastObj3 < thisObj  1.0e5 * CoinMax(fabs(thisObj), fabs(lastObj3))  1.0e7 1262 1263 #if ABC_NORMAL_DEBUG >01248 && firstFree_ < 0 && thisInf >= lastInf) { 1249 #if ABC_NORMAL_DEBUG > 0 1264 1250 if (handler_>logLevel() == 63) 1265 1251 printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_); 1266 1252 #endif 1267 1253 int maxFactor = abcFactorization_>maximumPivots(); 1268 1254 if (maxFactor > 10) { 1269 1270 1271 1272 #if ABC_NORMAL_DEBUG >01273 1274 1275 #endif 1276 } 1277 } else if (lastInf < testValue  trueInfeasibility == 1.123456e10) {1255 if (forceFactorization_ < 0) 1256 forceFactorization_ = maxFactor; 1257 forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3); 1258 #if ABC_NORMAL_DEBUG > 0 1259 if (handler_>logLevel() == 63) 1260 printf("Reducing factorization frequency to %d\n", forceFactorization_); 1261 #endif 1262 } 1263 } else if (lastInf < testValue  trueInfeasibility == 1.123456e10) { 1278 1264 if (infeasibilityCost_ < 1.0e14) { 1279 1280 1281 1282 #if ABC_NORMAL_DEBUG >01283 1284 1285 #endif 1286 1265 infeasibilityCost_ *= 1.5; 1266 // reset looping criterion 1267 abcProgress_.reset(); 1268 #if ABC_NORMAL_DEBUG > 0 1269 if (handler_>logLevel() == 63) 1270 printf("increasing weight to %g\n", infeasibilityCost_); 1271 #endif 1272 gutsOfPrimalSolution(3); 1287 1273 } 1288 1274 } … … 1292 1278 #if CLP_CAUTION 1293 1279 // If twice nearly there ... 1294 if (lastAverageInfeasibility<2.0*dualTolerance_) { 1295 double averageInfeasibility = sumDualInfeasibilities_ / 1296 static_cast<double>(numberDualInfeasibilities_ + 1); 1297 printf("last av %g now %g\n",lastAverageInfeasibility, 1298 averageInfeasibility); 1299 if (averageInfeasibility<2.0*dualTolerance_) 1280 if (lastAverageInfeasibility < 2.0 * dualTolerance_) { 1281 double averageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 1); 1282 printf("last av %g now %g\n", lastAverageInfeasibility, 1283 averageInfeasibility); 1284 if (averageInfeasibility < 2.0 * dualTolerance_) 1300 1285 sumOfRelaxedDualInfeasibilities_ = 0.0; 1301 1286 } 1302 1287 #endif 1303 1288 // give code benefit of doubt 1304 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && 1305 sumOfRelaxedPrimalInfeasibilities_ == 0.0) { 1289 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) { 1306 1290 // say been feasible 1307 1291 abcState_ = CLP_ABC_BEEN_FEASIBLE; … … 1318 1302 double error = CoinMin(1.0e2, largestPrimalError_); 1319 1303 // allow tolerance at least slightly bigger than standard 1320 relaxedToleranceP = relaxedToleranceP + 1304 relaxedToleranceP = relaxedToleranceP + error; 1321 1305 int ninfeas = abcNonLinearCost_>numberInfeasibilities(); 1322 1306 double sum = abcNonLinearCost_>sumInfeasibilities(); 1323 double average = sum / static_cast< double>(ninfeas);1324 #if ABC_NORMAL_DEBUG >31307 double average = sum / static_cast< double >(ninfeas); 1308 #if ABC_NORMAL_DEBUG > 3 1325 1309 if (handler_>logLevel() > 0) 1326 1327 1310 printf("nonLinearCost says infeasible %d summing to %g\n", 1311 ninfeas, sum); 1328 1312 #endif 1329 1313 if (average > relaxedToleranceP) { 1330 1331 1332 1333 #if ABC_NORMAL_DEBUG >31334 1335 #endif 1336 1337 1338 #if ABC_NORMAL_DEBUG >31339 1340 1314 sumOfRelaxedPrimalInfeasibilities_ = sum; 1315 numberPrimalInfeasibilities_ = ninfeas; 1316 sumPrimalInfeasibilities_ = sum; 1317 #if ABC_NORMAL_DEBUG > 3 1318 bool unflagged = 1319 #endif 1320 unflag(); 1321 abcProgress_.clearTimesFlagged(); 1322 #if ABC_NORMAL_DEBUG > 3 1323 if (unflagged && handler_>logLevel() > 0) 1324 printf("  but flagged variables\n"); 1341 1325 #endif 1342 1326 } … … 1347 1331 if ((dualFeasible()  problemStatus_ == 4) && !ifValuesPass) { 1348 1332 // see if extra helps 1349 if (abcNonLinearCost_>numberInfeasibilities() && 1350 (abcNonLinearCost_>sumInfeasibilities() > 1.0e3  sumOfRelaxedPrimalInfeasibilities_) 1351 && !alwaysOptimal) { 1333 if (abcNonLinearCost_>numberInfeasibilities() && (abcNonLinearCost_>sumInfeasibilities() > 1.0e3  sumOfRelaxedPrimalInfeasibilities_) 1334 && !alwaysOptimal) { 1352 1335 //may need infeasiblity cost changed 1353 1336 // we can see if we can construct a ray … … 1355 1338 // put nonbasics to bounds in case tolerance moved 1356 1339 // put back original costs 1357 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);; 1340 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1341 ; 1358 1342 abcNonLinearCost_>checkInfeasibilities(0.0); 1359 1343 gutsOfPrimalSolution(3); 1360 1344 // put back original costs 1361 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);; 1345 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1346 ; 1362 1347 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 1363 1348 gutsOfPrimalSolution(3); 1364 1349 // may have fixed infeasibilities  double check 1365 1350 if (abcNonLinearCost_>numberInfeasibilities() == 0) { 1366 1367 1368 1351 // carry on 1352 problemStatus_ = 1; 1353 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 1369 1354 } else { 1370 if ((infeasibilityCost_ >= 1.0e18  1371 numberDualInfeasibilities_ == 0) && perturbation_ == 101 1372 && (specialOptions_&8192)==0) { 1373 goToDual = unPerturb(); // stop any further perturbation 1355 if ((infeasibilityCost_ >= 1.0e18  numberDualInfeasibilities_ == 0) && perturbation_ == 101 1356 && (specialOptions_ & 8192) == 0) { 1357 goToDual = unPerturb(); // stop any further perturbation 1374 1358 #ifndef TRY_ABC_GUS 1375 1376 1377 #endif 1378 1379 1380 } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e2 1359 if (abcNonLinearCost_>sumInfeasibilities() > 1.0e1) 1360 goToDual = false; 1361 #endif 1362 numberDualInfeasibilities_ = 1; // carry on 1363 problemStatus_ = 1; 1364 } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e2 1381 1365 #ifndef TRY_ABC_GUS 1382 &&((moreSpecialOptions_ & 256) == 0&&(specialOptions_ & 8192) == 0)1366 && ((moreSpecialOptions_ & 256) == 0 && (specialOptions_ & 8192) == 0) 1383 1367 #else 1384 &&(specialOptions_ & 8192) == 0 1385 #endif 1386 ) { 1387 goToDual = true; 1388 abcFactorization_>pivotTolerance(CoinMax(0.5, abcFactorization_>pivotTolerance())); 1389 } 1390 if (!goToDual) { 1391 if (infeasibilityCost_ >= 1.0e20  1392 numberDualInfeasibilities_ == 0) { 1393 // we are infeasible  use as ray 1394 delete [] ray_; 1395 ray_ = new double [numberRows_]; 1396 CoinMemcpyN(dual_, numberRows_, ray_); 1397 // and get feasible duals 1398 infeasibilityCost_ = 0.0; 1399 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);; 1400 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 1401 gutsOfPrimalSolution(3); 1402 // so will exit 1403 infeasibilityCost_ = 1.0e30; 1404 // reset infeasibilities 1405 sumPrimalInfeasibilities_ = abcNonLinearCost_>sumInfeasibilities();; 1406 numberPrimalInfeasibilities_ = 1407 abcNonLinearCost_>numberInfeasibilities(); 1408 } 1409 if (infeasibilityCost_ < 1.0e20) { 1410 infeasibilityCost_ *= 5.0; 1411 // reset looping criterion 1412 abcProgress_.reset(); 1413 changeMade_++; // say change made 1414 handler_>message(CLP_PRIMAL_WEIGHT, messages_) 1415 << infeasibilityCost_ 1416 << CoinMessageEol; 1417 // put back original costs and then check 1418 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);; 1419 abcNonLinearCost_>checkInfeasibilities(0.0); 1420 gutsOfPrimalSolution(3); 1421 problemStatus_ = 1; //continue 1368 && (specialOptions_ & 8192) == 0 1369 #endif 1370 ) { 1371 goToDual = true; 1372 abcFactorization_>pivotTolerance(CoinMax(0.5, abcFactorization_>pivotTolerance())); 1373 } 1374 if (!goToDual) { 1375 if (infeasibilityCost_ >= 1.0e20  numberDualInfeasibilities_ == 0) { 1376 // we are infeasible  use as ray 1377 delete[] ray_; 1378 ray_ = new double[numberRows_]; 1379 CoinMemcpyN(dual_, numberRows_, ray_); 1380 // and get feasible duals 1381 infeasibilityCost_ = 0.0; 1382 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1383 ; 1384 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 1385 gutsOfPrimalSolution(3); 1386 // so will exit 1387 infeasibilityCost_ = 1.0e30; 1388 // reset infeasibilities 1389 sumPrimalInfeasibilities_ = abcNonLinearCost_>sumInfeasibilities(); 1390 ; 1391 numberPrimalInfeasibilities_ = abcNonLinearCost_>numberInfeasibilities(); 1392 } 1393 if (infeasibilityCost_ < 1.0e20) { 1394 infeasibilityCost_ *= 5.0; 1395 // reset looping criterion 1396 abcProgress_.reset(); 1397 changeMade_++; // say change made 1398 handler_>message(CLP_PRIMAL_WEIGHT, messages_) 1399 << infeasibilityCost_ 1400 << CoinMessageEol; 1401 // put back original costs and then check 1402 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1403 ; 1404 abcNonLinearCost_>checkInfeasibilities(0.0); 1405 gutsOfPrimalSolution(3); 1406 problemStatus_ = 1; //continue 1422 1407 #ifndef TRY_ABC_GUS 1423 1408 goToDual = false; 1424 1409 #else 1425 if((specialOptions_&8192)==0&&!sumOfRelaxedDualInfeasibilities_)1426 goToDual=true;1427 #endif 1428 1429 1430 1431 1432 1410 if ((specialOptions_ & 8192) == 0 && !sumOfRelaxedDualInfeasibilities_) 1411 goToDual = true; 1412 #endif 1413 } else { 1414 // say infeasible 1415 problemStatus_ = 1; 1416 } 1417 } 1433 1418 } 1434 1419 } else { 1435 1420 // may be optimal 1436 1421 if (perturbation_ == 101) { 1437 1422 goToDual = unPerturb(); // stop any further perturbation 1438 1423 #ifndef TRY_ABC_GUS 1439 1424 if ((numberRows_ > 20000  numberDualInfeasibilities_) && !numberTimesOptimal_) 1440 1425 #else 1441 if ((specialOptions_&8192)!=0)1442 #endif 1443 1444 1426 if ((specialOptions_ & 8192) != 0) 1427 #endif 1428 goToDual = false; // Better to carry on a bit longer 1429 lastCleaned_ = 1; // carry on 1445 1430 } 1446 1431 bool unflagged = (unflag() != 0); 1447 1432 abcProgress_.clearTimesFlagged(); 1448 if ( 1449 1450 1451 1452 1453 1454 1455 1456 1457 1458 1459 1460 1461 1462 1463 1464 1465 1433 if (lastCleaned_ != numberIterations_  unflagged) { 1434 handler_>message(CLP_PRIMAL_OPTIMAL, messages_) 1435 << primalTolerance_ 1436 << CoinMessageEol; 1437 if (numberTimesOptimal_ < 4) { 1438 numberTimesOptimal_++; 1439 changeMade_++; // say change made 1440 if (numberTimesOptimal_ == 1) { 1441 // better to have small tolerance even if slower 1442 abcFactorization_>zeroTolerance(CoinMin(abcFactorization_>zeroTolerance(), 1.0e15)); 1443 } 1444 lastCleaned_ = numberIterations_; 1445 if (primalTolerance_ != dblParam_[ClpPrimalTolerance]) 1446 handler_>message(CLP_PRIMAL_ORIGINAL, messages_) 1447 << CoinMessageEol; 1448 double oldTolerance = primalTolerance_; 1449 primalTolerance_ = dblParam_[ClpPrimalTolerance]; 1450 // put back original costs and then check 1466 1451 #if 0 //ndef NDEBUG 1467 1452 double largestDifference=0.0; … … 1471 1456 printf("largest change in cost %g\n",largestDifference); 1472 1457 #endif 1473 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);;1474 abcNonLinearCost_>checkInfeasibilities(oldTolerance);1475 gutsOfPrimalSolution(3);1476 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && 1477 1478 1479 1480 1481 1482 1483 goToDual=false;1484 1485 1486 1487 1488 1489 1490 1491 1492 1493 1494 1495 1458 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1459 ; 1460 abcNonLinearCost_>checkInfeasibilities(oldTolerance); 1461 gutsOfPrimalSolution(3); 1462 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) { 1463 // say optimal (with these bounds etc) 1464 numberDualInfeasibilities_ = 0; 1465 sumDualInfeasibilities_ = 0.0; 1466 numberPrimalInfeasibilities_ = 0; 1467 sumPrimalInfeasibilities_ = 0.0; 1468 goToDual = false; 1469 } 1470 if (dualFeasible() && !abcNonLinearCost_>numberInfeasibilities() && lastCleaned_ >= 0) 1471 problemStatus_ = 0; 1472 else 1473 problemStatus_ = 1; 1474 } else { 1475 problemStatus_ = 0; // optimal 1476 if (lastCleaned_ < numberIterations_) { 1477 handler_>message(CLP_SIMPLEX_GIVINGUP, messages_) 1478 << CoinMessageEol; 1479 } 1480 } 1496 1481 } else { 1497 1498 1499 1500 1482 if (!alwaysOptimal  !sumOfRelaxedPrimalInfeasibilities_) 1483 problemStatus_ = 0; // optimal 1484 else 1485 problemStatus_ = 1; // infeasible 1501 1486 } 1502 1487 } … … 1505 1490 if (problemStatus_ == 5) { 1506 1491 if (abcNonLinearCost_>numberInfeasibilities()) { 1507 if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) { 1508 // back off weight 1509 infeasibilityCost_ = 1.0e13; 1510 // reset looping criterion 1511 abcProgress_.reset(); 1512 unPerturb(); // stop any further perturbation 1513 } 1514 //we need infeasiblity cost changed 1515 if (infeasibilityCost_ < 1.0e20) { 1516 infeasibilityCost_ *= 5.0; 1517 // reset looping criterion 1518 abcProgress_.reset(); 1519 changeMade_++; // say change made 1520 handler_>message(CLP_PRIMAL_WEIGHT, messages_) 1521 << infeasibilityCost_ 1522 << CoinMessageEol; 1523 // put back original costs and then check 1524 CoinAbcMemcpy(abcCost_,costSaved_,numberTotal_);; 1525 gutsOfPrimalSolution(3); 1526 problemStatus_ = 1; //continue 1527 } else { 1528 // say infeasible 1529 problemStatus_ = 1; 1530 // we are infeasible  use as ray 1531 delete [] ray_; 1532 ray_ = new double [numberRows_]; 1533 CoinMemcpyN(dual_, numberRows_, ray_); 1534 } 1492 if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) { 1493 // back off weight 1494 infeasibilityCost_ = 1.0e13; 1495 // reset looping criterion 1496 abcProgress_.reset(); 1497 unPerturb(); // stop any further perturbation 1498 } 1499 //we need infeasiblity cost changed 1500 if (infeasibilityCost_ < 1.0e20) { 1501 infeasibilityCost_ *= 5.0; 1502 // reset looping criterion 1503 abcProgress_.reset(); 1504 changeMade_++; // say change made 1505 handler_>message(CLP_PRIMAL_WEIGHT, messages_) 1506 << infeasibilityCost_ 1507 << CoinMessageEol; 1508 // put back original costs and then check 1509 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_); 1510 ; 1511 gutsOfPrimalSolution(3); 1512 problemStatus_ = 1; //continue 1513 } else { 1514 // say infeasible 1515 problemStatus_ = 1; 1516 // we are infeasible  use as ray 1517 delete[] ray_; 1518 ray_ = new double[numberRows_]; 1519 CoinMemcpyN(dual_, numberRows_, ray_); 1520 } 1535 1521 } else { 1536 1537 1522 // say unbounded 1523 problemStatus_ = 2; 1538 1524 } 1539 1525 } else { 1540 1526 // carry on 1541 1527 problemStatus_ = 1; 1542 if(type == 3 && !ifValuesPass) { 1543 //bool unflagged = 1544 unflag(); 1545 abcProgress_.clearTimesFlagged(); 1546 if (sumDualInfeasibilities_ < 1.0e3  1547 (sumDualInfeasibilities_ / static_cast<double> (numberDualInfeasibilities_)) < 1.0e5  1548 abcProgress_.lastIterationNumber(0) == numberIterations_) { 1549 if (!numberPrimalInfeasibilities_) { 1550 if (numberTimesOptimal_ < 4) { 1551 numberTimesOptimal_++; 1552 changeMade_++; // say change made 1553 } else { 1554 problemStatus_ = 0; 1555 secondaryStatus_ = 5; 1556 } 1557 } 1558 } 1528 if (type == 3 && !ifValuesPass) { 1529 //bool unflagged = 1530 unflag(); 1531 abcProgress_.clearTimesFlagged(); 1532 if (sumDualInfeasibilities_ < 1.0e3  (sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_)) < 1.0e5  abcProgress_.lastIterationNumber(0) == numberIterations_) { 1533 if (!numberPrimalInfeasibilities_) { 1534 if (numberTimesOptimal_ < 4) { 1535 numberTimesOptimal_++; 1536 changeMade_++; // say change made 1537 } else { 1538 problemStatus_ = 0; 1539 secondaryStatus_ = 5; 1540 } 1541 } 1542 } 1559 1543 } 1560 1544 } … … 1562 1546 if (problemStatus_ == 0) { 1563 1547 double objVal = (abcNonLinearCost_>feasibleCost() 1564 1548 + objective_>nonlinearOffset()); 1565 1549 objVal /= (objectiveScale_ * rhsScale_); 1566 1550 double tol = 1.0e10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e8; 1567 1551 if (fabs(objVal  objectiveValue_) > tol) { 1568 #if ABC_NORMAL_DEBUG >31552 #if ABC_NORMAL_DEBUG > 3 1569 1553 if (handler_>logLevel() > 0) 1570 1571 1554 printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n", 1555 objVal, objectiveValue_); 1572 1556 #endif 1573 1557 objectiveValue_ = objVal; … … 1593 1577 moreSpecialOptions_ &= ~16; // clear check accuracy flag 1594 1578 #ifndef TRY_ABC_GUS 1595 if ((moreSpecialOptions_ & 256) != 0 saveSum(specialOptions_ & 8192) != 0)1596 goToDual =false;1579 if ((moreSpecialOptions_ & 256) != 0  saveSum  (specialOptions_ & 8192) != 0) 1580 goToDual = false; 1597 1581 #else 1598 1582 if ((specialOptions_ & 8192) != 0) 1599 goToDual=false; 1600 #endif 1601 if (goToDual  (numberIterations_ > 1000+baseIteration_ && largestPrimalError_ > 1.0e6 1602 && largestDualError_ > 1.0e6)) { 1583 goToDual = false; 1584 #endif 1585 if (goToDual  (numberIterations_ > 1000 + baseIteration_ && largestPrimalError_ > 1.0e6 && largestDualError_ > 1.0e6)) { 1603 1586 problemStatus_ = 10; // try dual 1604 1587 // See if second call 1605 if ((moreSpecialOptions_ & 256) != 0 abcNonLinearCost_>sumInfeasibilities()>1.0e20) {1588 if ((moreSpecialOptions_ & 256) != 0  abcNonLinearCost_>sumInfeasibilities() > 1.0e20) { 1606 1589 numberPrimalInfeasibilities_ = abcNonLinearCost_>numberInfeasibilities(); 1607 1590 sumPrimalInfeasibilities_ = abcNonLinearCost_>sumInfeasibilities(); 1608 1591 // say infeasible 1609 if (numberPrimalInfeasibilities_ &&(abcState_&CLP_ABC_BEEN_FEASIBLE)==0)1610 1592 if (numberPrimalInfeasibilities_ && (abcState_ & CLP_ABC_BEEN_FEASIBLE) == 0) 1593 problemStatus_ = 1; 1611 1594 } 1612 1595 } … … 1617 1600 } 1618 1601 #ifdef TRY_SPLIT_VALUES_PASS 1619 if (valuesChunk>0) { 1620 bool doFixing=firstFree_<0; 1621 if (numberIterations_==baseIteration_&& 1622 numberDualInfeasibilitiesWithoutFree_+1000<numberDualInfeasibilities_) { 1623 doFixing=true; 1624 int numberSuperBasic=0; 1602 if (valuesChunk > 0) { 1603 bool doFixing = firstFree_ < 0; 1604 if (numberIterations_ == baseIteration_ && numberDualInfeasibilitiesWithoutFree_ + 1000 < numberDualInfeasibilities_) { 1605 doFixing = true; 1606 int numberSuperBasic = 0; 1625 1607 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { 1626 1627 1628 } 1629 keepValuesPass =static_cast<int>((1.01*numberSuperBasic)/valuesChunk);1630 doOrdinaryVariables =static_cast<int>(valuesRatio*keepValuesPass);1631 } else if (valuesStop >0) {1632 if (numberIterations_ >=valuesStopproblemStatus_>=0) {1633 1634 1635 1636 int numberSuperBasic=0;1637 1638 1639 if (abcUpper_[iSequence]>abcLower_[iSequence]) {1640 setInternalStatus(iSequence,superBasic);1641 1642 1643 1644 1645 1646 1647 problemStatus_=1;1648 1649 1650 doFixing=false;1651 1652 valuesStop=1;1608 if (getInternalStatus(iSequence) == superBasic) 1609 numberSuperBasic++; 1610 } 1611 keepValuesPass = static_cast< int >((1.01 * numberSuperBasic) / valuesChunk); 1612 doOrdinaryVariables = static_cast< int >(valuesRatio * keepValuesPass); 1613 } else if (valuesStop > 0) { 1614 if (numberIterations_ >= valuesStop  problemStatus_ >= 0) { 1615 gutsOfSolution(3); 1616 abcNonLinearCost_>refreshFromPerturbed(primalTolerance_); 1617 gutsOfSolution(3); 1618 int numberSuperBasic = 0; 1619 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { 1620 if (getInternalStatus(iSequence) == isFixed) { 1621 if (abcUpper_[iSequence] > abcLower_[iSequence]) { 1622 setInternalStatus(iSequence, superBasic); 1623 numberSuperBasic++; 1624 } 1625 } 1626 } 1627 if (numberSuperBasic) { 1628 stateOfProblem_ = VALUES_PASS; 1629 problemStatus_ = 1; 1630 gutsOfSolution(3); 1631 } else { 1632 doFixing = false; 1633 } 1634 valuesStop = 1; 1653 1635 } else { 1654 doFixing=false;1655 } 1656 } 1636 doFixing = false; 1637 } 1638 } 1657 1639 if (doFixing) { 1658 1640 abcNonLinearCost_>refreshFromPerturbed(primalTolerance_); 1659 int numberSuperBasic=0;1641 int numberSuperBasic = 0; 1660 1642 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { 1661 1662 if (abcUpper_[iSequence]>abcLower_[iSequence]) 1663 setInternalStatus(iSequence,superBasic);1664 1665 if (getInternalStatus(iSequence) == superBasic) 1666 1667 } 1668 int numberFixed =0;1669 firstFree_ =1;1643 if (getInternalStatus(iSequence) == isFixed) { 1644 if (abcUpper_[iSequence] > abcLower_[iSequence]) 1645 setInternalStatus(iSequence, superBasic); 1646 } 1647 if (getInternalStatus(iSequence) == superBasic) 1648 numberSuperBasic++; 1649 } 1650 int numberFixed = 0; 1651 firstFree_ = 1; 1670 1652 if (numberSuperBasic) { 1671 double threshold = static_cast<double>(keepValuesPass)/numberSuperBasic;1672 if (1.1*keepValuesPass>numberSuperBasic) 1673 threshold=1.1;1674 1675 1676 if (randomNumberGenerator_.randomDouble()<threshold) {1677 if (firstFree_<0)1678 firstFree_=iSequence;1679 1680 1681 abcLower_[iSequence]=abcSolution_[iSequence];1682 abcUpper_[iSequence]=abcSolution_[iSequence];1683 setInternalStatus(iSequence,isFixed);1684 1685 1686 1653 double threshold = static_cast< double >(keepValuesPass) / numberSuperBasic; 1654 if (1.1 * keepValuesPass > numberSuperBasic) 1655 threshold = 1.1; 1656 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) { 1657 if (getInternalStatus(iSequence) == superBasic) { 1658 if (randomNumberGenerator_.randomDouble() < threshold) { 1659 if (firstFree_ < 0) 1660 firstFree_ = iSequence; 1661 } else { 1662 numberFixed++; 1663 abcLower_[iSequence] = abcSolution_[iSequence]; 1664 abcUpper_[iSequence] = abcSolution_[iSequence]; 1665 setInternalStatus(iSequence, isFixed); 1666 } 1667 } 1668 } 1687 1669 } 1688 1670 if (!numberFixed) { 1689 1690 valuesChunk=1;1691 } 1692 } 1693 } 1694 #endif 1695 if (doFactorization doWeights) {1671 // end 1672 valuesChunk = 1; 1673 } 1674 } 1675 } 1676 #endif 1677 if (doFactorization  doWeights) { 1696 1678 // restore weights (if saved)  also recompute infeasibility list 1697 1679 if (tentativeStatus > 3) … … 1702 1684 } 1703 1685 // Computes solutions  1 do duals, 2 do primals, 3 both 1704 int 1705 AbcSimplex::gutsOfPrimalSolution(int type) 1686 int AbcSimplex::gutsOfPrimalSolution(int type) 1706 1687 { 1707 1688 //work space 1708 1689 int whichArray[2]; 1709 for (int i =0;i<2;i++)1710 whichArray[i] =getAvailableArray();1711 CoinIndexedVector * 1712 CoinIndexedVector * 1690 for (int i = 0; i < 2; i++) 1691 whichArray[i] = getAvailableArray(); 1692 CoinIndexedVector *array1 = &usefulArray_[whichArray[0]]; 1693 CoinIndexedVector *array2 = &usefulArray_[whichArray[1]]; 1713 1694 // do work and check 1714 int numberRefinements =0;1715 if ((type &2)!=0) {1716 numberRefinements =computePrimals(array1,array2);1695 int numberRefinements = 0; 1696 if ((type & 2) != 0) { 1697 numberRefinements = computePrimals(array1, array2); 1717 1698 if (algorithm_ > 0 && abcNonLinearCost_ != NULL) { 1718 1699 // primal algorithm … … 1720 1701 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 1721 1702 if (abcNonLinearCost_>numberInfeasibilities()) 1722 1723 1724 1725 1726 1727 1703 if (handler_>detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) { 1704 handler_>message(CLP_SIMPLEX_NONLINEAR, messages_) 1705 << abcNonLinearCost_>changeInCost() 1706 << abcNonLinearCost_>numberInfeasibilities() 1707 << CoinMessageEol; 1708 } 1728 1709 } 1729 1710 checkPrimalSolution(true); 1730 1711 } 1731 if ((type &1)!=01732 #if CAN_HAVE_ZERO_OBJ >11733 &&(specialOptions_&2097152)==01734 #endif 1735 1736 numberRefinements +=computeDuals(NULL,array1,array2);1712 if ((type & 1) != 0 1713 #if CAN_HAVE_ZERO_OBJ > 1 1714 && (specialOptions_ & 2097152) == 0 1715 #endif 1716 ) { 1717 numberRefinements += computeDuals(NULL, array1, array2); 1737 1718 checkDualSolution(); 1738 1719 } 1739 for (int i =0;i<2;i++)1720 for (int i = 0; i < 2; i++) 1740 1721 setAvailableArray(whichArray[i]); 1741 rawObjectiveValue_ += sumNonBasicCosts_;1722 rawObjectiveValue_ += sumNonBasicCosts_; 1742 1723 //computeObjective(); // ? done in checkDualSolution?? 1743 1724 setClpSimplexObjectiveValue(); 1744 if (handler_>logLevel() > 3  (largestPrimalError_ > 1.0e2  1745 largestDualError_ > 1.0e2)) 1725 if (handler_>logLevel() > 3  (largestPrimalError_ > 1.0e2  largestDualError_ > 1.0e2)) 1746 1726 handler_>message(CLP_SIMPLEX_ACCURACY, messages_) 1747 1727 << largestPrimalError_ 1748 1728 << largestDualError_ 1749 1729 << CoinMessageEol; 1750 if ((largestPrimalError_ > 1.0e1 largestDualError_>1.0e1)1751 1752 #if ABC_NORMAL_DEBUG >01730 if ((largestPrimalError_ > 1.0e1  largestDualError_ > 1.0e1) 1731 && numberRows_ > 100 && numberIterations_) { 1732 #if ABC_NORMAL_DEBUG > 0 1753 1733 printf("Large errors  primal %g dual %g\n", 1754 largestPrimalError_,largestDualError_);1734 largestPrimalError_, largestDualError_); 1755 1735 #endif 1756 1736 // Change factorization tolerance … … 1767 1747 On exit rhsArray will have changes in costs of basic variables 1768 1748 */ 1769 void 1770 AbcSimplexPrimal::primalRow(CoinIndexedVector * rowArray, 1771 CoinIndexedVector * rhsArray, 1772 CoinIndexedVector * spareArray, 1773 int valuesPass) 1749 void AbcSimplexPrimal::primalRow(CoinIndexedVector *rowArray, 1750 CoinIndexedVector *rhsArray, 1751 CoinIndexedVector *spareArray, 1752 int valuesPass) 1774 1753 { 1775 1754 #if 1 1776 for (int iRow =0;iRow<numberRows_;iRow++) {1777 int iPivot =abcPivotVariable_[iRow];1778 assert (costBasic_[iRow]==abcCost_[iPivot]);1779 assert (lowerBasic_[iRow]==abcLower_[iPivot]);1780 assert (upperBasic_[iRow]==abcUpper_[iPivot]);1781 assert (solutionBasic_[iRow]==abcSolution_[iPivot]);1755 for (int iRow = 0; iRow < numberRows_; iRow++) { 1756 int iPivot = abcPivotVariable_[iRow]; 1757 assert(costBasic_[iRow] == abcCost_[iPivot]); 1758 assert(lowerBasic_[iRow] == abcLower_[iPivot]); 1759 assert(upperBasic_[iRow] == abcUpper_[iPivot]); 1760 assert(solutionBasic_[iRow] == abcSolution_[iPivot]); 1782 1761 } 1783 1762 #endif … … 1785 1764 if (valuesPass) { 1786 1765 dualIn_ = abcCost_[sequenceIn_]; 1787 1788 double * 1766 1767 double *work = rowArray>denseVector(); 1789 1768 int number = rowArray>getNumElements(); 1790 int * 1791 1769 int *which = rowArray>getIndices(); 1770 1792 1771 int iIndex; 1793 1772 for (iIndex = 0; iIndex < number; iIndex++) { 1794 1773 1795 1774 int iRow = which[iIndex]; 1796 1775 double alpha = work[iRow]; … … 1805 1784 // towards nearest bound 1806 1785 if (valueIn_  lowerIn_ < upperIn_  valueIn_) { 1807 1808 1786 directionIn_ = 1; 1787 dualIn_ = dualTolerance_; 1809 1788 } else { 1810 1811 1812 } 1813 } 1814 } 1815 1789 directionIn_ = 1; 1790 dualIn_ = dualTolerance_; 1791 } 1792 } 1793 } 1794 1816 1795 // sequence stays as row number until end 1817 1796 pivotRow_ = 1; 1818 1797 int numberRemaining = 0; 1819 1798 1820 1799 double totalThru = 0.0; // for when variables flip 1821 1800 // Allow first few iterations to take tiny … … 1833 1812 double lastPivot = 0.0; 1834 1813 double lastTheta = 1.0e50; 1835 1814 1836 1815 // use spareArrays to put ones looked at in 1837 1816 // First one is list of candidates … … 1839 1818 // Second array has current set of pivot candidates 1840 1819 // with a backup list saved in double * part of indexed vector 1841 1820 1842 1821 // pivot elements 1843 double * 1822 double *spare; 1844 1823 // indices 1845 int * 1824 int *index; 1846 1825 spareArray>clear(); 1847 1826 spare = spareArray>denseVector(); 1848 1827 index = spareArray>getIndices(); 1849 1828 1850 1829 // we also need somewhere for effective rhs 1851 double * 1830 double *rhs = rhsArray>denseVector(); 1852 1831 // and we can use indices to point to alpha 1853 1832 // that way we can store fabs(alpha) 1854 int * 1833 int *indexPoint = rhsArray>getIndices(); 1855 1834 //int numberFlip=0; // Those which may change if flips 1856 1835 1857 1836 /* 1858 1837 First we get a list of possible pivots. We can also see if the … … 1864 1843 1865 1844 */ 1866 1845 1867 1846 // do first pass to get possibles 1868 1847 // We can also see if unbounded 1869 1870 double * 1848 1849 double *work = rowArray>denseVector(); 1871 1850 int number = rowArray>getNumElements(); 1872 int * 1873 1851 int *which = rowArray>getIndices(); 1852 1874 1853 // we need to swap sign if coming in from ub 1875 1854 double way = directionIn_; … … 1879 1858 else 1880 1859 maximumMovement = CoinMin(1.0e30, valueIn_  lowerIn_); 1881 1860 1882 1861 double averageTheta = abcNonLinearCost_>averageTheta(); 1883 1862 double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement); … … 1892 1871 // but make a bit more pessimistic 1893 1872 dualCheck = CoinMax(dualCheck  100.0 * dualTolerance_, 0.99 * dualCheck); 1894 1873 1895 1874 int iIndex; 1896 1875 //#define CLP_DEBUG 1897 #if ABC_NORMAL_DEBUG >31876 #if ABC_NORMAL_DEBUG > 3 1898 1877 if (numberIterations_ == 1369  numberIterations_ == 3840) { 1899 1878 double dj = abcCost_[sequenceIn_]; 1900 1879 printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_); 1901 1880 for (iIndex = 0; iIndex < number; iIndex++) { 1902 1881 1903 1882 int iRow = which[iIndex]; 1904 1883 double alpha = work[iRow]; … … 1906 1885 dj = alpha * costBasic_[iRow]; 1907 1886 printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n", 1908 1909 1887 iRow, iPivot, lowerBasic_[iRow], solutionBasic_[iRow], upperBasic_[iRow], 1888 alpha, solutionBasic_[iRow]  1.0e9 * alpha, costBasic_[iRow], dj); 1910 1889 } 1911 1890 } … … 1920 1899 #endif 1921 1900 for (iIndex = 0; iIndex < number; iIndex++) { 1922 1901 1923 1902 int iRow = which[iIndex]; 1924 1903 double alpha = work[iRow]; … … 1931 1910 // do computation same way as later on in primal 1932 1911 if (alpha > 0.0) { 1933 1934 1935 1936 1937 1938 1912 // basic variable going towards lower bound 1913 double bound = lowerBasic_[iRow]; 1914 // must be exactly same as when used 1915 double change = tentativeTheta * alpha; 1916 possible = (oldValue  change) <= bound + primalTolerance_; 1917 oldValue = bound; 1939 1918 } else { 1940 1941 1942 1943 1944 1945 1946 alpha = alpha;1919 // basic variable going towards upper bound 1920 double bound = upperBasic_[iRow]; 1921 // must be exactly same as when used 1922 double change = tentativeTheta * alpha; 1923 possible = (oldValue  change) >= bound  primalTolerance_; 1924 oldValue = bound  oldValue; 1925 alpha = alpha; 1947 1926 } 1948 1927 double value; 1949 1928 //assert (oldValue >= 10.0*tolerance); 1950 1929 if (possible) { 1951 1930 value = oldValue  upperTheta * alpha; 1952 1931 #ifdef CLP_USER_DRIVEN1 1953 if(!userChoiceValid1(this,iPivot,oldValue,1954 upperTheta,alpha,work[iRow]*way))1955 value =0.0; // say can't use1956 #endif 1957 1958 1959 1960 1961 1962 1963 1964 1965 1966 1967 1968 1969 1970 } 1971 } 1972 if (upperTheta < maximumMovement && totalThru *infeasibilityCost_ >= 1.0001 * dualCheck) {1932 if (!userChoiceValid1(this, iPivot, oldValue, 1933 upperTheta, alpha, work[iRow] * way)) 1934 value = 0.0; // say can't use 1935 #endif 1936 if (value < primalTolerance_ && alpha >= acceptablePivot) { 1937 upperTheta = (oldValue + primalTolerance_) / alpha; 1938 } 1939 // add to list 1940 spare[numberRemaining] = alpha; 1941 rhs[numberRemaining] = oldValue; 1942 indexPoint[numberRemaining] = iIndex; 1943 index[numberRemaining++] = iRow; 1944 totalThru += alpha; 1945 setActive(iRow); 1946 //} else if (value<primalTolerance_*1.002) { 1947 // May change if is a flip 1948 //indexRhs[numberFlip++]=iRow; 1949 } 1950 } 1951 if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) { 1973 1952 // Can pivot here 1974 1953 break; … … 1982 1961 } 1983 1962 totalThru = 0.0; 1984 1963 1985 1964 theta_ = maximumMovement; 1986 1965 1987 1966 bool goBackOne = false; 1988 1967 if (objective_>type() > 1) 1989 1968 dualIn_ = saveDj; 1990 1969 1991 1970 //printf("%d remain out of %d\n",numberRemaining,number); 1992 1971 int iTry = 0; … … 1994 1973 if (numberRemaining && upperTheta < maximumMovement) { 1995 1974 // First check if previously chosen one will work 1996 1975 1997 1976 // first get ratio with tolerance 1998 for ( 1999 1977 for (; iTry < MAXTRY; iTry++) { 1978 2000 1979 upperTheta = maximumMovement; 2001 1980 int iBest = 1; 2002 1981 for (iIndex = 0; iIndex < numberRemaining; iIndex++) { 2003 2004 2005 2006 2007 1982 1983 double alpha = spare[iIndex]; 1984 double oldValue = rhs[iIndex]; 1985 double value = oldValue  upperTheta * alpha; 1986 2008 1987 #ifdef CLP_USER_DRIVEN1 2009 int sequenceOut=abcPivotVariable_[index[iIndex]];2010 if(!userChoiceValid1(this,sequenceOut,oldValue,2011 upperTheta,alpha, 0.0))2012 value =0.0; // say can't use2013 #endif 2014 2015 2016 2017 2018 } 2019 1988 int sequenceOut = abcPivotVariable_[index[iIndex]]; 1989 if (!userChoiceValid1(this, sequenceOut, oldValue, 1990 upperTheta, alpha, 0.0)) 1991 value = 0.0; // say can't use 1992 #endif 1993 if (value < primalTolerance_ && alpha >= acceptablePivot) { 1994 upperTheta = (oldValue + primalTolerance_) / alpha; 1995 iBest = iIndex; // just in case weird numbers 1996 } 1997 } 1998 2020 1999 // now look at best in this lot 2021 2000 // But also see how infeasible small pivots will make … … 2024 2003 pivotRow_ = 1; 2025 2004 for (iIndex = 0; iIndex < numberRemaining; iIndex++) { 2026 2027 2028 2029 2030 2031 2032 2033 2034 2035 2036 2037 2038 2039 2040 2041 2005 2006 int iRow = index[iIndex]; 2007 double alpha = spare[iIndex]; 2008 double oldValue = rhs[iIndex]; 2009 double value = oldValue  upperTheta * alpha; 2010 2011 if (value <= 0  iBest == iIndex) { 2012 // how much would it cost to go thru and modify bound 2013 double trueAlpha = way * work[iRow]; 2014 totalThru += abcNonLinearCost_>changeInCost(iRow, trueAlpha, rhs[iIndex]); 2015 setActive(iRow); 2016 if (alpha > bestPivot) { 2017 bestPivot = alpha; 2018 theta_ = oldValue / bestPivot; 2019 pivotRow_ = iRow; 2020 } else if (alpha < acceptablePivot 2042 2021 #ifdef CLP_USER_DRIVEN1 2043 !userChoiceValid1(this,abcPivotVariable_[index[iIndex]], 2044 oldValue,upperTheta,alpha,0.0) 2045 #endif 2046 ) { 2047 if (value < primalTolerance_) 2048 sumInfeasibilities += value  primalTolerance_; 2049 } 2050 } 2051 } 2052 if (bestPivot < 0.1 * bestEverPivot && 2053 bestEverPivot > 1.0e6 && bestPivot < 1.0e3) { 2054 // back to previous one 2055 goBackOne = true; 2056 break; 2022  !userChoiceValid1(this, abcPivotVariable_[index[iIndex]], 2023 oldValue, upperTheta, alpha, 0.0) 2024 #endif 2025 ) { 2026 if (value < primalTolerance_) 2027 sumInfeasibilities += value  primalTolerance_; 2028 } 2029 } 2030 } 2031 if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e6 && bestPivot < 1.0e3) { 2032 // back to previous one 2033 goBackOne = true; 2034 break; 2057 2035 } else if (pivotRow_ == 1 && upperTheta > largeValue_) { 2058 2059 2060 2061 2062 2063 2064 2036 if (lastPivot > acceptablePivot) { 2037 // back to previous one 2038 goBackOne = true; 2039 } else { 2040 // can only get here if all pivots so far too small 2041 } 2042 break; 2065 2043 } else if (totalThru >= dualCheck) { 2066 2067 2068 2069 2070 2071 2072 2073 2074 2075 2044 if (sumInfeasibilities > primalTolerance_ && !abcNonLinearCost_>numberInfeasibilities()) { 2045 // Looks a bad choice 2046 if (lastPivot > acceptablePivot) { 2047 goBackOne = true; 2048 } else { 2049 // say no good 2050 dualIn_ = 0.0; 2051 } 2052 } 2053 break; // no point trying another loop 2076 2054 } else { 2077 2078 2079 2080 2055 lastPivotRow = pivotRow_; 2056 lastTheta = theta_; 2057 if (bestPivot > bestEverPivot) 2058 bestEverPivot = bestPivot; 2081 2059 } 2082 2060 } … … 2095 2073 tentativeTheta = 1.0e50; 2096 2074 for (iIndex = 0; iIndex < number; iIndex++) { 2097 2098 2099 2100 2101 2102 2103 2104 2105 2106 2107 2108 2109 2110 2111 alpha = alpha;2112 2113 2114 2115 2075 int iRow = which[iIndex]; 2076 double alpha = work[iRow]; 2077 alpha *= way; 2078 double oldValue = solutionBasic_[iRow]; 2079 // get where in bound sequence 2080 // note that after this alpha is actually fabs(alpha) 2081 if (alpha > 0.0) { 2082 // basic variable going towards lower bound 2083 double bound = lowerBasic_[iRow]; 2084 oldValue = bound; 2085 } else { 2086 // basic variable going towards upper bound 2087 double bound = upperBasic_[iRow]; 2088 oldValue = bound  oldValue; 2089 alpha = alpha; 2090 } 2091 if (oldValue  tentativeTheta * alpha < 0.0) { 2092 tentativeTheta = oldValue / alpha; 2093 } 2116 2094 } 2117 2095 // If free in then see if we can get to 0.0 2118 2096 if (lowerIn_ < 1.0e20 && upperIn_ > 1.0e20) { 2119 2120 2121 2122 2123 2097 if (dualIn_ * valueIn_ > 0.0) { 2098 if (fabs(valueIn_) < 1.0e2 && (tentativeTheta < fabs(valueIn_)  tentativeTheta > 1.0e20)) { 2099 tentativeTheta = fabs(valueIn_); 2100 } 2101 } 2124 2102 } 2125 2103 if (tentativeTheta < 1.0e10) 2126 2104 valueOut_ = valueIn_ + way * tentativeTheta; 2127 2105 } 2128 2106 } … … 2158 2136 theta_ = minimumTheta; 2159 2137 for (iIndex = 0; iIndex < numberRemaining  numberRemaining; iIndex++) { 2160 2161 2138 largestInfeasibility = CoinMax(largestInfeasibility, 2139 (rhs[iIndex]  spare[iIndex] * theta_)); 2162 2140 } 2163 2141 //#define CLP_DEBUG 2164 #if ABC_NORMAL_DEBUG >32142 #if ABC_NORMAL_DEBUG > 3 2165 2143 if (largestInfeasibility > primalTolerance_ && (handler_>logLevel() & 32) > 1) 2166 2167 2144 printf("Primal tolerance increased from %g to %g\n", 2145 primalTolerance_, largestInfeasibility); 2168 2146 #endif 2169 2147 //#undef CLP_DEBUG … … 2173 2151 if (theta_ > tentativeTheta) { 2174 2152 for (iIndex = 0; iIndex < number; iIndex++) 2175 2153 setActive(which[iIndex]); 2176 2154 } 2177 2155 if (way < 0.0) 2178 theta_ =  2156 theta_ = theta_; 2179 2157 double newValue = valueOut_  theta_ * alpha_; 2180 2158 // If 4 bit set  Force outgoing variables to exact bound (primal) 2181 2159 if (alpha_ * way < 0.0) { 2182 directionOut_ = 1; 2160 directionOut_ = 1; // to upper bound 2183 2161 if (fabs(theta_) > 1.0e6  (specialOptions_ & 4) != 0) { 2184 2162 upperOut_ = abcNonLinearCost_>nearest(pivotRow_, newValue); 2185 2163 } else { 2186 2164 upperOut_ = newValue; 2187 2165 } 2188 2166 } else { 2189 directionOut_ = 1; 2167 directionOut_ = 1; // to lower bound 2190 2168 if (fabs(theta_) > 1.0e6  (specialOptions_ & 4) != 0) { 2191 2169 lowerOut_ = abcNonLinearCost_>nearest(pivotRow_, newValue); 2192 2170 } else { 2193 2171 lowerOut_ = newValue; 2194 2172 } 2195 2173 } … … 2205 2183 alpha_ = 0.0; 2206 2184 if (way < 0.0) { 2207 directionOut_ = 1; 2185 directionOut_ = 1; // to lower bound 2208 2186 theta_ = lowerOut_  valueOut_; 2209 2187 } else { 2210 directionOut_ = 1; 2188 directionOut_ = 1; // to upper bound 2211 2189 theta_ = upperOut_  valueOut_; 2212 2190 } 2213 2191 } 2214 2192 2215 2193 double theta1 = CoinMax(theta_, 1.0e12); 2216 2194 double theta2 = numberIterations_ * abcNonLinearCost_>averageTheta(); 2217 2195 // Set average theta 2218 abcNonLinearCost_>setAverageTheta((theta1 + theta2) / (static_cast< double>(numberIterations_ + 1)));2196 abcNonLinearCost_>setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1))); 2219 2197 //if (numberIterations_%1000==0) 2220 2198 //printf("average theta is %g\n",abcNonLinearCost_>averageTheta()); 2221 2199 2222 2200 // clear arrays 2223 2201 CoinZeroN(spare, numberRemaining); 2224 2202 CoinZeroN(rhs, numberRemaining); 2225 2203 2226 2204 // put back original bounds etc 2227 2205 CoinMemcpyN(index, numberRemaining, rhsArray>getIndices()); … … 2230 2208 abcNonLinearCost_>goBackAll(rhsArray); 2231 2209 rhsArray>clear(); 2232 2233 2210 } 2234 2211 /* … … 2239 2216 On exit rhsArray will have changes in costs of basic variables 2240 2217 */ 2241 void 2242 AbcSimplexPrimal::primalRow(CoinIndexedVector * rowArray, 2243 CoinIndexedVector * rhsArray, 2244 CoinIndexedVector * spareArray, 2245 pivotStruct & stuff) 2218 void AbcSimplexPrimal::primalRow(CoinIndexedVector *rowArray, 2219 CoinIndexedVector *rhsArray, 2220 CoinIndexedVector *spareArray, 2221 pivotStruct &stuff) 2246 2222 { 2247 int valuesPass =stuff.valuesPass_;2248 double dualIn =stuff.dualIn_;2249 double lowerIn =stuff.lowerIn_;2250 double upperIn =stuff.upperIn_;2251 double valueIn =stuff.valueIn_;2252 int sequenceIn =stuff.sequenceIn_;2253 int directionIn =stuff.directionIn_;2254 int pivotRow =1;2255 int sequenceOut =1;2256 double dualOut =0.0;2257 double lowerOut =0.0;2258 double upperOut =0.0;2259 double valueOut =0.0;2223 int valuesPass = stuff.valuesPass_; 2224 double dualIn = stuff.dualIn_; 2225 double lowerIn = stuff.lowerIn_; 2226 double upperIn = stuff.upperIn_; 2227 double valueIn = stuff.valueIn_; 2228 int sequenceIn = stuff.sequenceIn_; 2229 int directionIn = stuff.directionIn_; 2230 int pivotRow = 1; 2231 int sequenceOut = 1; 2232 double dualOut = 0.0; 2233 double lowerOut = 0.0; 2234 double upperOut = 0.0; 2235 double valueOut = 0.0; 2260 2236 double theta; 2261 2237 double alpha; 2262 int directionOut =0;2238 int directionOut = 0; 2263 2239 if (valuesPass) { 2264 2240 dualIn = abcCost_[sequenceIn]; 2265 2266 double * 2241 2242 double *work = rowArray>denseVector(); 2267 2243 int number = rowArray>getNumElements(); 2268 int * 2269 2244 int *which = rowArray>getIndices(); 2245 2270 2246 int iIndex; 2271 2247 for (iIndex = 0; iIndex < number; iIndex++) { 2272 2248 2273 2249 int iRow = which[iIndex]; 2274 2250 double alpha = work[iRow]; … … 2283 2259 // towards nearest bound 2284 2260 if (valueIn  lowerIn < upperIn  valueIn) { 2285 2286 2261 directionIn = 1; 2262 dualIn = dualTolerance_; 2287 2263 } else { 2288 2289 2290 } 2291 } 2292 } 2293 2264 directionIn = 1; 2265 dualIn = dualTolerance_; 2266 } 2267 } 2268 } 2269 2294 2270 // sequence stays as row number until end 2295 2271 pivotRow = 1; 2296 2272 int numberRemaining = 0; 2297 2273 2298 2274 double totalThru = 0.0; // for when variables flip 2299 2275 // Allow first few iterations to take tiny … … 2311 2287 double lastPivot = 0.0; 2312 2288 double lastTheta = 1.0e50; 2313 2289 2314 2290 // use spareArrays to put ones looked at in 2315 2291 // First one is list of candidates … … 2317 2293 // Second array has current set of pivot candidates 2318 2294 // with a backup list saved in double * part of indexed vector 2319 2295 2320 2296 // pivot elements 2321 double * 2297 double *spare; 2322 2298 // indices 2323 int * 2299 int *index; 2324 2300 spareArray>clear(); 2325 2301 spare = spareArray>denseVector(); 2326 2302 index = spareArray>getIndices(); 2327 2303 2328 2304 // we also need somewhere for effective rhs 2329 double * 2305 double *rhs = rhsArray>denseVector(); 2330 2306 // and we can use indices to point to alpha 2331 2307 // that way we can store fabs(alpha) 2332 int * 2308 int *indexPoint = rhsArray>getIndices(); 2333 2309 //int numberFlip=0; // Those which may change if flips 2334 2310 2335 2311 /* 2336 2312 First we get a list of possible pivots. We can also see if the … … 2342 2318 2343 2319 */ 2344 2320 2345 2321 // do first pass to get possibles 2346 2322 // We can also see if unbounded 2347 2348 double * 2323 2324 double *work = rowArray>denseVector(); 2349 2325 int number = rowArray>getNumElements(); 2350 int * 2351 2326 int *which = rowArray>getIndices(); 2327 2352 2328 // we need to swap sign if coming in from ub 2353 2329 double way = directionIn; … … 2357 2333 else 2358 2334 maximumMovement = CoinMin(1.0e30, valueIn  lowerIn); 2359 2335 2360 2336 double averageTheta = abcNonLinearCost_>averageTheta(); 2361 2337 double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement); … … 2370 2346 // but make a bit more pessimistic 2371 2347 dualCheck = CoinMax(dualCheck  100.0 * dualTolerance_, 0.99 * dualCheck); 2372 2348 2373 2349 int iIndex; 2374 2350 while (true) { … … 2388 2364 // do computation same way as later on in primal 2389 2365 if (alpha > 0.0) { 2390 2391 2392 2393 2394 2395 2366 // basic variable going towards lower bound 2367 double bound = lowerBasic_[iRow]; 2368 // must be exactly same as when used 2369 double change = tentativeTheta * alpha; 2370 possible = (oldValue  change) <= bound + primalTolerance_; 2371 oldValue = bound; 2396 2372 } else { 2397 2398 2399 2400 2401 2402 2403 alpha = alpha;2373 // basic variable going towards upper bound 2374 double bound = upperBasic_[iRow]; 2375 // must be exactly same as when used 2376 double change = tentativeTheta * alpha; 2377 possible = (oldValue  change) >= bound  primalTolerance_; 2378 oldValue = bound  oldValue; 2379 alpha = alpha; 2404 2380 } 2405 2381 double value; 2406 2382 if (possible) { 2407 2383 value = oldValue  upperTheta * alpha; 2408 2384 #ifdef CLP_USER_DRIVEN1 2409 if(!userChoiceValid1(this,iPivot,oldValue,2410 upperTheta,alpha,work[iRow]*way))2411 value =0.0; // say can't use2412 #endif 2413 2414 2415 2416 2417 2418 2419 2420 2421 2422 2423 } 2424 } 2425 if (upperTheta < maximumMovement && totalThru *infeasibilityCost_ >= 1.0001 * dualCheck) {2385 if (!userChoiceValid1(this, iPivot, oldValue, 2386 upperTheta, alpha, work[iRow] * way)) 2387 value = 0.0; // say can't use 2388 #endif 2389 if (value < primalTolerance_ && alpha >= acceptablePivot) { 2390 upperTheta = (oldValue + primalTolerance_) / alpha; 2391 } 2392 // add to list 2393 spare[numberRemaining] = alpha; 2394 rhs[numberRemaining] = oldValue; 2395 indexPoint[numberRemaining] = iIndex; 2396 index[numberRemaining++] = iRow; 2397 totalThru += alpha; 2398 setActive(iRow); 2399 } 2400 } 2401 if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) { 2426 2402 // Can pivot here 2427 2403 break; … … 2435 2411 } 2436 2412 totalThru = 0.0; 2437 2413 2438 2414 theta = maximumMovement; 2439 2415 2440 2416 bool goBackOne = false; 2441 2417 2442 2418 //printf("%d remain out of %d\n",numberRemaining,number); 2443 2419 int iTry = 0; 2444 2420 if (numberRemaining && upperTheta < maximumMovement) { 2445 2421 // First check if previously chosen one will work 2446 2422 2447 2423 // first get ratio with tolerance 2448 for ( 2449 2424 for (; iTry < MAXTRY; iTry++) { 2425 2450 2426 upperTheta = maximumMovement; 2451 2427 int iBest = 1; 2452 2428 for (iIndex = 0; iIndex < numberRemaining; iIndex++) { 2453 2454 2455 2456 2457 2429 2430 double alpha = spare[iIndex]; 2431 double oldValue = rhs[iIndex]; 2432 double value = oldValue  upperTheta * alpha; 2433 2458 2434 #ifdef CLP_USER_DRIVEN1 2459 int sequenceOut=abcPivotVariable_[index[iIndex]];2460 if(!userChoiceValid1(this,sequenceOut,oldValue,2461 upperTheta,alpha, 0.0))2462 value =0.0; // say can't use2463 #endif 2464 2465 2466 2467 2468 } 2469 2435 int sequenceOut = abcPivotVariable_[index[iIndex]]; 2436 if (!userChoiceValid1(this, sequenceOut, oldValue, 2437 upperTheta, alpha, 0.0)) 2438 value = 0.0; // say can't use 2439 #endif 2440 if (value < primalTolerance_ && alpha >= acceptablePivot) { 2441 upperTheta = (oldValue + primalTolerance_) / alpha; 2442 iBest = iIndex; // just in case weird numbers 2443 } 2444 } 2445 2470 2446 // now look at best in this lot 2471 2447 // But also see how infeasible small pivots will make … … 2474 2450 pivotRow = 1; 2475 2451 for (iIndex = 0; iIndex < numberRemaining; iIndex++) { 2476 2477 2478 2479 2480 2481 2482 2483 2484 2485 2486 2487 2488 2489 2490 2491 2452 2453 int iRow = index[iIndex]; 2454 double alpha = spare[iIndex]; 2455 double oldValue = rhs[iIndex]; 2456 double value = oldValue  upperTheta * alpha; 2457 2458 if (value <= 0  iBest == iIndex) { 2459 // how much would it cost to go thru and modify bound 2460 double trueAlpha = way * work[iRow]; 2461 totalThru += abcNonLinearCost_>changeInCost(iRow, trueAlpha, rhs[iIndex]); 2462 setActive(iRow); 2463 if (alpha > bestPivot) { 2464 bestPivot = alpha; 2465 theta = oldValue / bestPivot; 2466 pivotRow = iRow; 2467 } else if (alpha < acceptablePivot 2492 2468 #ifdef CLP_USER_DRIVEN1 2493 !userChoiceValid1(this,abcPivotVariable_[index[iIndex]], 2494 oldValue,upperTheta,alpha,0.0) 2495 #endif 2496 ) { 2497 if (value < primalTolerance_) 2498 sumInfeasibilities += value  primalTolerance_; 2499 } 2500 } 2501 } 2502 if (bestPivot < 0.1 * bestEverPivot && 2503 bestEverPivot > 1.0e6 && bestPivot < 1.0e3) { 2504 // back to previous one 2505 goBackOne = true; 2506 break; 2469  !userChoiceValid1(this, abcPivotVariable_[index[iIndex]], 2470 oldValue, upperTheta, alpha, 0.0) 2471 #endif 2472 ) { 2473 if (value < primalTolerance_) 2474 sumInfeasibilities += value  primalTolerance_; 2475 } 2476 } 2477 } 2478 if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e6 && bestPivot < 1.0e3) { 2479 // back to previous one 2480 goBackOne = true; 2481 break; 2507 2482 } else if (pivotRow == 1 && upperTheta > largeValue_) { 2508 2509 2510 2511 2512 2513 2514 2483 if (lastPivot > acceptablePivot) { 2484 // back to previous one 2485 goBackOne = true; 2486 } else { 2487 // can only get here if all pivots so far too small 2488 } 2489 break; 2515 2490 } else if (totalThru >= dualCheck) { 2516 2517 2518 2519 2520 2521 2522 2523 2524 2525 2491 if (sumInfeasibilities > primalTolerance_ && !abcNonLinearCost_>numberInfeasibilities()) { 2492 // Looks a bad choice 2493 if (lastPivot > acceptablePivot) { 2494 goBackOne = true; 2495 } else { 2496 // say no good 2497 dualIn = 0.0; 2498 } 2499 } 2500 break; // no point trying another loop 2526 2501 } else { 2527 2528 2529 2530 2502 lastPivotRow = pivotRow; 2503 lastTheta = theta; 2504 if (bestPivot > bestEverPivot) 2505 bestEverPivot = bestPivot; 2531 2506 } 2532 2507 } … … 2545 2520 tentativeTheta = 1.0e50; 2546 2521 for (iIndex = 0; iIndex < number; iIndex++) { 2547 2548 2549 2550 2551 2552 2553 2554 2555 2556 2557 2558 2559 2560 2561 alpha = alpha;2562 2563 2564 2565 2522 int iRow = which[iIndex]; 2523 double alpha = work[iRow]; 2524 alpha *= way; 2525 double oldValue = solutionBasic_[iRow]; 2526 // get where in bound sequence 2527 // note that after this alpha is actually fabs(alpha) 2528 if (alpha > 0.0) { 2529 // basic variable going towards lower bound 2530 double bound = lowerBasic_[iRow]; 2531 oldValue = bound; 2532 } else { 2533 // basic variable going towards upper bound 2534 double bound = upperBasic_[iRow]; 2535 oldValue = bound  oldValue; 2536 alpha = alpha; 2537 } 2538 if (oldValue  tentativeTheta * alpha < 0.0) { 2539 tentativeTheta = oldValue / alpha; 2540 } 2566 2541 } 2567 2542 // If free in then see if we can get to 0.0 2568 2543 if (lowerIn < 1.0e20 && upperIn > 1.0e20) { 2569 2570 2571 2572 2573 2544 if (dualIn * valueIn > 0.0) { 2545 if (fabs(valueIn) < 1.0e2 && (tentativeTheta < fabs(valueIn)  tentativeTheta > 1.0e20)) { 2546 tentativeTheta = fabs(valueIn); 2547 } 2548 } 2574 2549 } 2575 2550 if (tentativeTheta < 1.0e10) 2576 2551 valueOut = valueIn + way * tentativeTheta; 2577 2552 } 2578 2553 } … … 2612 2587 if (theta > tentativeTheta) { 2613 2588 for (iIndex = 0; iIndex < number; iIndex++) 2614 2589 setActive(which[iIndex]); 2615 2590 } 2616 2591 if (way < 0.0) 2617 theta =  2592 theta = theta; 2618 2593 double newValue = valueOut  theta * alpha; 2619 2594 // If 4 bit set  Force outgoing variables to exact bound (primal) 2620 2595 if (alpha * way < 0.0) { 2621 directionOut = 1; 2596 directionOut = 1; // to upper bound 2622 2597 if (fabs(theta) > 1.0e6  (specialOptions_ & 4) != 0) { 2623 2598 upperOut = abcNonLinearCost_>nearest(pivotRow, newValue); 2624 2599 } else { 2625 2600 upperOut = newValue; 2626 2601 } 2627 2602 } else { 2628 directionOut = 1; 2603 directionOut = 1; // to lower bound 2629 2604 if (fabs(theta) > 1.0e6  (specialOptions_ & 4) != 0) { 2630 2605 lowerOut = abcNonLinearCost_>nearest(pivotRow, newValue); 2631 2606 } else { 2632 2607 lowerOut = newValue; 2633 2608 } 2634 2609 } … … 2644 2619 alpha = 0.0; 2645 2620 if (way < 0.0) { 2646 directionOut = 1; 2621 directionOut = 1; // to lower bound 2647 2622 theta = lowerOut  valueOut; 2648 2623 } else { 2649 directionOut = 1; 2624 directionOut = 1; // to upper bound 2650 2625 theta = upperOut  valueOut; 2651 2626 } 2652 2627 } 2653 2654 2628 2655 2629 // clear arrays 2656 2630 CoinZeroN(spare, numberRemaining); 2657 2631 CoinZeroN(rhs, numberRemaining); 2658 2632 2659 2633 // put back original bounds etc 2660 2634 CoinMemcpyN(index, numberRemaining, rhsArray>getIndices()); … … 2662 2636 abcNonLinearCost_>goBackAll(rhsArray); 2663 2637 rhsArray>clear(); 2664 stuff.theta_ =theta;2665 stuff.alpha_ =alpha;2666 stuff.dualIn_ =dualIn;2667 stuff.dualOut_ =dualOut;2668 stuff.lowerOut_ =lowerOut;2669 stuff.upperOut_ =upperOut;2670 stuff.valueOut_ =valueOut;2671 stuff.sequenceOut_ =sequenceOut;2672 stuff.directionOut_ =directionOut;2673 stuff.pivotRow_ =pivotRow;2638 stuff.theta_ = theta; 2639 stuff.alpha_ = alpha; 2640 stuff.dualIn_ = dualIn; 2641 stuff.dualOut_ = dualOut; 2642 stuff.lowerOut_ = lowerOut; 2643 stuff.upperOut_ = upperOut; 2644 stuff.valueOut_ = valueOut; 2645 stuff.sequenceOut_ = sequenceOut; 2646 stuff.directionOut_ = directionOut; 2647 stuff.pivotRow_ = pivotRow; 2674 2648 } 2675 2649 /* … … 2680 2654 For easy problems we can just choose one of the first columns we look at 2681 2655 */ 2682 void 2683 AbcSimplexPrimal::primalColumn(CoinPartitionedVector * updates, 2684 CoinPartitionedVector * spareRow2, 2685 CoinPartitionedVector * spareColumn1) 2656 void AbcSimplexPrimal::primalColumn(CoinPartitionedVector *updates, 2657 CoinPartitionedVector *spareRow2, 2658 CoinPartitionedVector *spareColumn1) 2686 2659 { 2687 for (int i =0;i<4;i++)2688 multipleSequenceIn_[i] =1;2660 for (int i = 0; i < 4; i++) 2661 multipleSequenceIn_[i] = 1; 2689 2662 sequenceIn_ = abcPrimalColumnPivot_>pivotColumn(updates, 2690 2663 spareRow2, spareColumn1); 2691 2664 if (sequenceIn_ >= 0) { 2692 2665 valueIn_ = abcSolution_[sequenceIn_]; … … 2706 2679 After rowArray will have list of cost changes 2707 2680 */ 2708 int 2709 AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray, 2710 double theta, 2711 double & objectiveChange, 2712 int valuesPass) 2681 int AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector *rowArray, 2682 double theta, 2683 double &objectiveChange, 2684 int valuesPass) 2713 2685 { 2714 2686 // Cost on pivot row may change  may need to change dualIn … … 2716 2688 if (pivotRow_ >= 0) 2717 2689 oldCost = abcCost_[sequenceOut_]; 2718 double * 2690 double *work = rowArray>denseVector(); 2719 2691 int number = rowArray>getNumElements(); 2720 int * 2721 2692 int *which = rowArray>getIndices(); 2693 2722 2694 int newNumber = 0; 2723 2695 abcNonLinearCost_>setChangeInCost(0.0); … … 2728 2700 if (!valuesPass) { 2729 2701 for (iIndex = 0; iIndex < number; iIndex++) { 2730 2702 2731 2703 int iRow = which[iIndex]; 2732 2704 double alpha = work[iRow]; … … 2741 2713 // check if not active then okay 2742 2714 if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != 1) { 2743 2744 2745 2746 2747 2748 2749 2750 2751 2752 2753 2754 2755 2756 2757 2758 2759 2760 2715 // But make sure one going out is feasible 2716 if (change > 0.0) { 2717 // going down 2718 if (value <= abcLower_[iPivot] + primalTolerance_) { 2719 if (iPivot == sequenceOut_ && value > abcLower_[iPivot]  relaxedTolerance) 2720 value = abcLower_[iPivot]; 2721 //double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2722 //assert (!difference  fabs(change) > 1.0e9); 2723 } 2724 } else { 2725 // going up 2726 if (value >= abcUpper_[iPivot]  primalTolerance_) { 2727 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance) 2728 value = abcUpper_[iPivot]; 2729 //double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2730 //assert (!difference  fabs(change) > 1.0e9); 2731 } 2732 } 2761 2733 } 2762 2734 #endif 2763 2735 if (active(iRow)  theta_ < 0.0) { 2764 2765 2766 2767 2768 2769 2770 2771 2772 2773 2774 2775 2776 2777 2778 2779 2780 2781 2782 2783 2784 2785 2786 2787 2788 2789 2790 2791 2792 2736 clearActive(iRow); 2737 // But make sure one going out is feasible 2738 if (change > 0.0) { 2739 // going down 2740 if (value <= abcLower_[iPivot] + primalTolerance_) { 2741 if (iPivot == sequenceOut_ && value >= abcLower_[iPivot]  relaxedTolerance) 2742 value = abcLower_[iPivot]; 2743 double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2744 if (difference) { 2745 work[iRow] = difference; 2746 //change reduced cost on this 2747 //abcDj_[iPivot] = difference; 2748 which[newNumber++] = iRow; 2749 } 2750 } 2751 } else { 2752 // going up 2753 if (value >= abcUpper_[iPivot]  primalTolerance_) { 2754 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance) 2755 value = abcUpper_[iPivot]; 2756 double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2757 if (difference) { 2758 work[iRow] = difference; 2759 //change reduced cost on this 2760 //abcDj_[iPivot] = difference; 2761 which[newNumber++] = iRow; 2762 } 2763 } 2764 } 2793 2765 } 2794 2766 } … … 2796 2768 // values pass so look at all 2797 2769 for (iIndex = 0; iIndex < number; iIndex++) { 2798 2770 2799 2771 int iRow = which[iIndex]; 2800 2772 double alpha = work[iRow]; … … 2809 2781 // But make sure one going out is feasible 2810 2782 if (change > 0.0) { 2811 2812 2813 2814 2815 2816 2817 2818 2819 2820 2821 2822 2783 // going down 2784 if (value <= abcLower_[iPivot] + primalTolerance_) { 2785 if (iPivot == sequenceOut_ && value > abcLower_[iPivot]  relaxedTolerance) 2786 value = abcLower_[iPivot]; 2787 double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2788 if (difference) { 2789 work[iRow] = difference; 2790 //change reduced cost on this 2791 abcDj_[iPivot] = difference; 2792 which[newNumber++] = iRow; 2793 } 2794 } 2823 2795 } else { 2824 2825 2826 2827 2828 2829 2830 2831 2832 2833 2834 2835 2796 // going up 2797 if (value >= abcUpper_[iPivot]  primalTolerance_) { 2798 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance) 2799 value = abcUpper_[iPivot]; 2800 double difference = abcNonLinearCost_>setOneBasic(iRow, value); 2801 if (difference) { 2802 work[iRow] = difference; 2803 //change reduced cost on this 2804 abcDj_[iPivot] = difference; 2805 which[newNumber++] = iRow; 2806 } 2807 } 2836 2808 } 2837 2809 } … … 2852 2824 } 2853 2825 // Perturbs problem 2854 void 2855 AbcSimplexPrimal::perturb(int /*type*/) 2826 void AbcSimplexPrimal::perturb(int /*type*/) 2856 2827 { 2857 2828 if (perturbation_ > 100) … … 2868 2839 double largestPositive; 2869 2840 matrix_>rangeOfElements(smallestNegative, largestNegative, 2870 2841 smallestPositive, largestPositive); 2871 2842 smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive); 2872 2843 largestPositive = CoinMax(fabs(largestNegative), largestPositive); 2873 2844 double elementRatio = largestPositive / smallestPositive; 2874 if ((!numberIterations_  initialSumInfeasibilities_==1.23456789e5)&& perturbation_ == 50) {2845 if ((!numberIterations_  initialSumInfeasibilities_ == 1.23456789e5) && perturbation_ == 50) { 2875 2846 // See if we need to perturb 2876 2847 int numberTotal = CoinMax(numberRows_, numberColumns_); 2877 double * 2848 double *sort = new double[numberTotal]; 2878 2849 int nFixed = 0; 2879 2850 for (i = 0; i < numberRows_; i++) { … … 2882 2853 double value = 0.0; 2883 2854 if (lo && lo < 1.0e20) { 2884 2885 2886 2887 2888 2889 2890 2855 if (up && up < 1.0e20) { 2856 value = 0.5 * (lo + up); 2857 if (lo == up) 2858 nFixed++; 2859 } else { 2860 value = lo; 2861 } 2891 2862 } else { 2892 2893 2863 if (up && up < 1.0e20) 2864 value = up; 2894 2865 } 2895 2866 sort[i] = value; … … 2900 2871 for (i = 1; i < numberRows_; i++) { 2901 2872 if (last != sort[i]) 2902 2873 number++; 2903 2874 last = sort[i]; 2904 2875 } … … 2909 2880 if (number * 4 > numberRows_  elementRatio > 1.0e12) { 2910 2881 perturbation_ = 100; 2911 delete 2882 delete[] sort; 2912 2883 // Make sure feasible bounds 2913 2884 if (abcNonLinearCost_) { 2914 2915 2916 2885 abcNonLinearCost_>refresh(); 2886 abcNonLinearCost_>checkInfeasibilities(); 2887 //abcNonLinearCost_>feasibleBounds(); 2917 2888 } 2918 2889 moveToBasic(); … … 2925 2896 nFixed = 0; 2926 2897 for (i = 0; i < numberColumns_; i++) { 2927 2928 2929 2930 2931 2932 2933 2934 2935 2936 2937 2938 2939 2940 2941 2942 2898 double lo = fabs(columnLower_[i]); 2899 double up = fabs(columnUpper_[i]); 2900 double value = 0.0; 2901 if (lo && lo < 1.0e20) { 2902 if (up && up < 1.0e20) { 2903 value = 0.5 * (lo + up); 2904 if (lo == up) 2905 nFixed++; 2906 } else { 2907 value = lo; 2908 } 2909 } else { 2910 if (up && up < 1.0e20) 2911 value = up; 2912 } 2913 sort[i] = value; 2943 2914 } 2944 2915 std::sort(sort, sort + numberColumns_); … … 2946 2917 last = sort[0]; 2947 2918 for (i = 1; i < numberColumns_; i++) { 2948 2949 2950 2919 if (last != sort[i]) 2920 number++; 2921 last = sort[i]; 2951 2922 } 2952 2923 //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_), … … 2954 2925 } 2955 2926 #endif 2956 delete 2927 delete[] sort; 2957 2928 if (number * 4 > numberColumns_) { 2958 2929 perturbation_ = 100; 2959 2930 // Make sure feasible bounds 2960 2931 if (abcNonLinearCost_) { 2961 2962 2963 2932 abcNonLinearCost_>refresh(); 2933 abcNonLinearCost_>checkInfeasibilities(); 2934 //abcNonLinearCost_>feasibleBounds(); 2964 2935 } 2965 2936 moveToBasic(); … … 2973 2944 // maximum fraction of rhs/bounds to perturb 2974 2945 double maximumFraction = 1.0e5; 2975 double overallMultiplier = (perturbation_==50perturbation_>54) ? 2.0 : 0.2;2946 double overallMultiplier = (perturbation_ == 50  perturbation_ > 54) ? 2.0 : 0.2; 2976 2947 #ifdef HEAVY_PERTURBATION 2977 if (perturbation_==50)2978 perturbation_=HEAVY_PERTURBATION;2948 if (perturbation_ == 50) 2949 perturbation_ = HEAVY_PERTURBATION; 2979 2950 #endif 2980 2951 if (perturbation_ >= 50) { … … 2982 2953 for (i = 0; i < numberColumns_ + numberRows_; i++) { 2983 2954 if (abcUpper_[i] > abcLower_[i] + primalTolerance_) { 2984 2985 2986 2987 2988 2989 2990 2991 2992 2993 2994 2955 double lowerValue, upperValue; 2956 if (abcLower_[i] > 1.0e20) 2957 lowerValue = fabs(abcLower_[i]); 2958 else 2959 lowerValue = 0.0; 2960 if (abcUpper_[i] < 1.0e20) 2961 upperValue = fabs(abcUpper_[i]); 2962 else 2963 upperValue = 0.0; 2964 double value = CoinMax(fabs(lowerValue), fabs(upperValue)); 2965 value = CoinMin(value, abcUpper_[i]  abcLower_[i]); 2995 2966 #if 1 2996 2997 2998 2999 2967 if (value) { 2968 perturbation += value; 2969 numberNonZero++; 2970 } 3000 2971 #else 3001 2972 perturbation = CoinMax(perturbation, value); 3002 2973 #endif 3003 2974 } 3004 2975 } 3005 2976 if (numberNonZero) 3006 perturbation /= static_cast< double>(numberNonZero);2977 perturbation /= static_cast< double >(numberNonZero); 3007 2978 else 3008 2979 perturbation = 1.0e1; … … 3010 2981 // reduce 3011 2982 while (perturbation_ < 55) { 3012 3013 3014 2983 perturbation_++; 2984 perturbation *= 0.25; 2985 bias *= 0.25; 3015 2986 } 3016 2987 perturbation_ = 50; … … 3018 2989 // increase 3019 2990 while (perturbation_ > 55) { 3020 3021 3022 2991 overallMultiplier *= 1.2; 2992 perturbation_; 2993 perturbation *= 4.0; 3023 2994 } 3024 2995 perturbation_ = 50; … … 3045 3016 maximumFraction *= 0.1; 3046 3017 } 3047 if (savePerturbation ==51) {3048 perturbation = CoinMin(0.1, perturbation);3049 maximumFraction *= 0.1;3018 if (savePerturbation == 51) { 3019 perturbation = CoinMin(0.1, perturbation); 3020 maximumFraction *= 0.1; 3050 3021 } 3051 3022 //if (number != numberRows_) … … 3058 3029 //const double * COIN_RESTRICT perturbationArray = perturbationSaved_; 3059 3030 // Make sure feasible bounds 3060 if (abcNonLinearCost_ &&true) {3031 if (abcNonLinearCost_ && true) { 3061 3032 abcNonLinearCost_>refresh(); 3062 3033 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); … … 3064 3035 } 3065 3036 double tolerance = 100.0 * primalTolerance_; 3066 int numberChanged =0;3037 int numberChanged = 0; 3067 3038 // Set bit if fixed 3068 for (int i =0;i<numberRows_;i++) {3069 if (rowLower_[i] !=rowUpper_[i])3039 for (int i = 0; i < numberRows_; i++) { 3040 if (rowLower_[i] != rowUpper_[i]) 3070 3041 internalStatus_[i] &= ~128; 3071 3042 else 3072 3043 internalStatus_[i] = 128; 3073 3044 } 3074 for (int i =0;i<numberColumns_;i++) {3075 if (columnLower_[i] !=columnUpper_[i])3076 internalStatus_[i +numberRows_] &= ~128;3045 for (int i = 0; i < numberColumns_; i++) { 3046 if (columnLower_[i] != columnUpper_[i]) 3047 internalStatus_[i + numberRows_] &= ~128; 3077 3048 else 3078 internalStatus_[i +numberRows_] = 128;3049 internalStatus_[i + numberRows_] = 128; 3079 3050 } 3080 3051 //double multiplier = perturbation*maximumFraction; 3081 3052 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) { 3082 3053 if (getInternalStatus(iSequence) == basic) { 3083 if ((internalStatus_[i] & 128)!=0)3084 3054 if ((internalStatus_[i] & 128) != 0) 3055 continue; 3085 3056 double lowerValue = abcLower_[iSequence]; 3086 3057 double upperValue = abcUpper_[iSequence]; 3087 3058 if (upperValue > lowerValue + tolerance) { 3088 3089 3090 3091 3092 3093 3094 value = CoinMax(value,primalTolerance_);3095 double perturbationValue=overallMultiplier*randomNumberGenerator_.randomDouble();3096 3097 3098 3099 3100 if (abcNonLinearCost_>getCurrentStatus(iSequence)==CLP_FEASIBLE)3101 abcPerturbation_[iSequence]=abcLower_[iSequence];3102 3103 abcPerturbation_[iSequence+numberTotal_]=abcLower_[iSequence];3104 3105 3106 3107 if (abcNonLinearCost_>getCurrentStatus(iSequence)==CLP_FEASIBLE)3108 abcPerturbation_[iSequence+numberTotal_]=abcUpper_[iSequence];3109 3110 abcPerturbation_[iSequence]=abcUpper_[iSequence];3111 3059 double solutionValue = abcSolution_[iSequence]; 3060 double difference = upperValue  lowerValue; 3061 difference = CoinMin(difference, perturbation); 3062 difference = CoinMin(difference, fabs(solutionValue) + 1.0); 3063 double value = maximumFraction * (difference + bias); 3064 value = CoinMin(value, 0.1); 3065 value = CoinMax(value, primalTolerance_); 3066 double perturbationValue = overallMultiplier * randomNumberGenerator_.randomDouble(); 3067 value *= perturbationValue; 3068 if (solutionValue  lowerValue <= primalTolerance_) { 3069 abcLower_[iSequence] = value; 3070 // change correct saved value 3071 if (abcNonLinearCost_>getCurrentStatus(iSequence) == CLP_FEASIBLE) 3072 abcPerturbation_[iSequence] = abcLower_[iSequence]; 3073 else 3074 abcPerturbation_[iSequence + numberTotal_] = abcLower_[iSequence]; 3075 } else if (upperValue  solutionValue <= primalTolerance_) { 3076 abcUpper_[iSequence] += value; 3077 // change correct saved value 3078 if (abcNonLinearCost_>getCurrentStatus(iSequence) == CLP_FEASIBLE) 3079 abcPerturbation_[iSequence + numberTotal_] = abcUpper_[iSequence]; 3080 else 3081 abcPerturbation_[iSequence] = abcUpper_[iSequence]; 3082 } else { 3112 3083 #if 0 3113 3084 if (iSequence >= numberColumns_) { … … 3121 3092 } 3122 3093 #else 3123 3124 #endif 3125 3126 3127 3128 3129 3130 3131 3132 3133 if (value > (fabs(solutionValue) + 1.0)*largestPerCent)3134 3135 3136 3137 3138 3094 value = 0.0; 3095 #endif 3096 } 3097 if (value) { 3098 numberChanged++; 3099 if (printOut) 3100 printf("col %d lower from %g to %g, upper from %g to %g\n", 3101 iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue); 3102 if (solutionValue) { 3103 largest = CoinMax(largest, value); 3104 if (value > (fabs(solutionValue) + 1.0) * largestPerCent) 3105 largestPerCent = value / (fabs(solutionValue) + 1.0); 3106 } else { 3107 largestZero = CoinMax(largestZero, value); 3108 } 3109 } 3139 3110 } 3140 3111 } … … 3144 3115 for (iSequence = 0; iSequence < maximumAbcNumberRows_ + numberColumns_; iSequence++) { 3145 3116 if (getInternalStatus(iSequence) != basic) { 3146 if ((internalStatus_[i] &128)!=0)3147 3148 3149 3150 3151 3152 3153 3154 3155 3156 3157 value = CoinMax(value,primalTolerance_);3158 double perturbationValue=overallMultiplier*randomNumberGenerator_.randomDouble();3159 3160 3161 3162 3163 if (abcNonLinearCost_>getCurrentStatus(iSequence)==CLP_FEASIBLE)3164 abcPerturbation_[iSequence]=abcLower_[iSequence];3165 3166 abcPerturbation_[iSequence+numberTotal_]=abcLower_[iSequence];3167 3168 3169 3170 if (abcNonLinearCost_>getCurrentStatus(iSequence)==CLP_FEASIBLE)3171 abcPerturbation_[iSequence+numberTotal_]=abcUpper_[iSequence];3172 3173 abcPerturbation_[iSequence]=abcUpper_[iSequence];3174 3175 3176 3177 3178 3179 3180 3181 3182 3183 if (value > (fabs(solutionValue) + 1.0)*largestPerCent)3184 3185 3186 3187 3188 3189 3117 if ((internalStatus_[i] & 128) != 0) 3118 continue; 3119 double lowerValue = abcLower_[iSequence]; 3120 double upperValue = abcUpper_[iSequence]; 3121 if (upperValue > lowerValue + tolerance) { 3122 double solutionValue = abcSolution_[iSequence]; 3123 double difference = upperValue  lowerValue; 3124 difference = CoinMin(difference, perturbation); 3125 difference = CoinMin(difference, fabs(solutionValue) + 1.0); 3126 double value = maximumFraction * (difference + bias); 3127 value = CoinMin(value, 0.1); 3128 value = CoinMax(value, primalTolerance_); 3129 double perturbationValue = overallMultiplier * randomNumberGenerator_.randomDouble(); 3130 value *= perturbationValue; 3131 if (solutionValue  lowerValue <= primalTolerance_) { 3132 abcLower_[iSequence] = value; 3133 // change correct saved value 3134 if (abcNonLinearCost_>getCurrentStatus(iSequence) == CLP_FEASIBLE) 3135 abcPerturbation_[iSequence] = abcLower_[iSequence]; 3136 else 3137 abcPerturbation_[iSequence + numberTotal_] = abcLower_[iSequence]; 3138 } else if (upperValue  solutionValue <= primalTolerance_) { 3139 abcUpper_[iSequence] += value; 3140 // change correct saved value 3141 if (abcNonLinearCost_>getCurrentStatus(iSequence) == CLP_FEASIBLE) 3142 abcPerturbation_[iSequence + numberTotal_] = abcUpper_[iSequence]; 3143 else 3144 abcPerturbation_[iSequence] = abcUpper_[iSequence]; 3145 } else { 3146 value = 0.0; 3147 } 3148 if (value) { 3149 if (printOut) 3150 printf("col %d lower from %g to %g, upper from %g to %g\n", 3151 iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue); 3152 if (solutionValue) { 3153 largest = CoinMax(largest, value); 3154 if (value > (fabs(solutionValue) + 1.0) * largestPerCent) 3155 largestPerCent = value / (fabs(solutionValue) + 1.0); 3156 } else { 3157 largestZero = CoinMax(largestZero, value); 3158 } 3159 } 3160 } 3190 3161 } 3191 3162 } … … 3194 3165 for (i = 0; i < numberColumns_ + numberRows_; i++) { 3195 3166 internalStatus_[i] &= ~128; 3196 switch (getInternalStatus(i)) {3197 3167 switch (getInternalStatus(i)) { 3168 3198 3169 case basic: 3199 3170 break; … … 3232 3203 } 3233 3204 // un perturb 3234 bool 3235 AbcSimplexPrimal::unPerturb() 3205 bool AbcSimplexPrimal::unPerturb() 3236 3206 { 3237 3207 if (perturbation_ != 101) … … 3240 3210 copyFromSaved(); 3241 3211 // copy bounds to perturbation 3242 CoinAbcMemcpy(abcPerturbation_, abcLower_,numberTotal_);3243 CoinAbcMemcpy(abcPerturbation_ +numberTotal_,abcUpper_,numberTotal_);3212 CoinAbcMemcpy(abcPerturbation_, abcLower_, numberTotal_); 3213 CoinAbcMemcpy(abcPerturbation_ + numberTotal_, abcUpper_, numberTotal_); 3244 3214 //sanityCheck(); 3245 3215 // unflag … … 3255 3225 abcNonLinearCost_>checkInfeasibilities(primalTolerance_); 3256 3226 // Try using dual 3257 return abcNonLinearCost_>sumInfeasibilities()!=0.0; 3258 3227 return abcNonLinearCost_>sumInfeasibilities() != 0.0; 3259 3228 } 3260 3229 // Unflag all variables and return number unflagged 3261 int 3262 AbcSimplexPrimal::unflag() 3230 int AbcSimplexPrimal::unflag() 3263 3231 { 3264 3232 int i; … … 3273 3241 // only say if reasonable dj 3274 3242 if (fabs(abcDj_[i]) > relaxedToleranceD) 3275 3276 } 3277 } 3278 #if ABC_NORMAL_DEBUG >03243 numberFlagged++; 3244 } 3245 } 3246 #if ABC_NORMAL_DEBUG > 0 3279 3247 if (handler_>logLevel() > 2 && numberFlagged && objective_>type() > 1) 3280 3248 printf("%d unflagged\n", numberFlagged); … … 3283 3251 } 3284 3252 // Do not change infeasibility cost and always say optimal 3285 void 3286 AbcSimplexPrimal::alwaysOptimal(bool onOff) 3253 void AbcSimplexPrimal::alwaysOptimal(bool onOff) 3287 3254 { 3288 3255 if (onOff) … … 3291 3258 specialOptions_ &= ~1; 3292 3259 } 3293 bool 3294 AbcSimplexPrimal::alwaysOptimal() const 3260 bool AbcSimplexPrimal::alwaysOptimal() const 3295 3261 { 3296 3262 return (specialOptions_ & 1) != 0; 3297 3263 } 3298 3264 // Flatten outgoing variables i.e.  always to exact bound 3299 void 3300 AbcSimplexPrimal::exactOutgoing(bool onOff) 3265 void AbcSimplexPrimal::exactOutgoing(bool onOff) 3301 3266 { 3302 3267 if (onOff) … … 3305 3270 specialOptions_ &= ~4; 3306 3271 } 3307 bool 3308 AbcSimplexPrimal::exactOutgoing() const 3272 bool AbcSimplexPrimal::exactOutgoing() const 3309 3273 { 3310 3274 return (specialOptions_ & 4) != 0; … … 3320 3284 +3 max iterations (iteration done) 3321 3285 */ 3322 int 3323 AbcSimplexPrimal::pivotResult(int ifValuesPass) 3286 int AbcSimplexPrimal::pivotResult(int ifValuesPass) 3324 3287 { 3325 3288 3326 3289 bool roundAgain = true; 3327 3290 int returnCode = 1; 3328 3291 3329 3292 // loop round if user setting and doing refactorization 3330 3293 while (roundAgain) { … … 3334 3297 sequenceOut_ = 1; 3335 3298 usefulArray_[arrayForFtran_].clear(); 3336 3299 3337 3300 // we found a pivot column 3338 3301 // update the incoming column … … 3345 3308 if ((moreSpecialOptions_ & 512) == 0) { 3346 3309 #endif 3347 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_], 3348 3349 3310 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_], 3311 &usefulArray_[arrayForTableauRow_], 3312 ifValuesPass); 3350 3313 #ifdef CLP_USER_DRIVEN 3351 3314 // user can tell which use it is 3352 3315 int status = eventHandler_>event(ClpEventHandler::pivotRow); 3353 3316 if (status >= 0) { 3354 3355 3356 3317 problemStatus_ = 5; 3318 secondaryStatus_ = ClpEventHandler::pivotRow; 3319 break; 3357 3320 } 3358 3321 } else { 3359 3322 int status = eventHandler_>event(ClpEventHandler::pivotRow); 3360 3323 if (status >= 0) { 3361 3362 3363 3324 problemStatus_ = 5; 3325 secondaryStatus_ = ClpEventHandler::pivotRow; 3326 break; 3364 3327 } 3365 3328 } … … 3369 3332 //assert (fabs(alpha_)>=1.0e5(objective_>type()<2!objective_>activated())pivotRow_==2); 3370 3333 if (pivotRow_ == 1  (pivotRow_ >= 0 && fabs(alpha_) < 1.0e5)) { 3371 if(fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_>type() < 2) {3372 3373 3374 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_], 3375 3376 3377 3378 3379 3380 3381 3382 3383 3384 3385 3386 3387 3388 3389 3390 3391 3392 3334 if (fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_>type() < 2) { 3335 // try other way 3336 directionIn_ = directionIn_; 3337 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_], 3338 &usefulArray_[arrayForTableauRow_], 3339 0); 3340 } 3341 if (pivotRow_ == 1  (pivotRow_ >= 0 && fabs(alpha_) < 1.0e5)) { 3342 // reject it 3343 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3344 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3345 << x << sequenceWithin(sequenceIn_) 3346 << CoinMessageEol; 3347 setFlagged(sequenceIn_); 3348 abcProgress_.incrementTimesFlagged(); 3349 abcProgress_.clearBadTimes(); 3350 lastBadIteration_ = numberIterations_; // say be more cautious 3351 clearAll(); 3352 pivotRow_ = 1; 3353 returnCode = 5; 3354 break; 3355 } 3393 3356 } 3394 3357 } … … 3406 3369 } 3407 3370 #endif 3408 if (!ifValuesPass && (saveDj * dualIn_ < test1  3409 fabs(saveDj  dualIn_) > checkValue*(1.0 + fabs(saveDj))  3410 fabs(dualIn_) < test2)) { 3411 if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) > 3412 1.0e5)) { 3413 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3414 handler_>message(CLP_PRIMAL_DJ, messages_) 3415 << x << sequenceWithin(sequenceIn_) 3416 << saveDj << dualIn_ 3417 << CoinMessageEol; 3418 if(lastGoodIteration_ != numberIterations_) { 3419 clearAll(); 3420 pivotRow_ = 1; // say no weights update 3421 returnCode = 4; 3422 if(lastGoodIteration_ + 1 == numberIterations_) { 3423 // not looking wonderful  try cleaning bounds 3424 // put nonbasics to bounds in case tolerance moved 3425 abcNonLinearCost_>checkInfeasibilities(0.0); 3426 } 3427 sequenceOut_ = 1; 3428 sequenceIn_=1; 3429 if (abcFactorization_>pivots()<10&&abcFactorization_>pivotTolerance()<0.25) 3430 abcFactorization_>saferTolerances(1.0,1.03); 3431 break; 3432 } else { 3433 // take on more relaxed criterion 3434 if (saveDj * dualIn_ < test1  3435 fabs(saveDj  dualIn_) > 2.0e1 * (1.0 + fabs(dualIn_))  3436 fabs(dualIn_) < test2) { 3437 // need to reject something 3438 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3439 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3440 << x << sequenceWithin(sequenceIn_) 3441 << CoinMessageEol; 3442 setFlagged(sequenceIn_); 3443 abcProgress_.incrementTimesFlagged(); 3444 #if 1 //def FEB_TRY 3445 // Make safer? 3446 double tolerance=abcFactorization_>pivotTolerance(); 3447 abcFactorization_>saferTolerances (1.0, 1.03); 3448 #endif 3449 abcProgress_.clearBadTimes(); 3450 lastBadIteration_ = numberIterations_; // say be more cautious 3451 clearAll(); 3452 pivotRow_ = 1; 3453 returnCode = 5; 3454 if(tolerance<abcFactorization_>pivotTolerance()) 3455 returnCode=4; 3456 sequenceOut_ = 1; 3457 sequenceIn_=1; 3458 break; 3459 } 3460 } 3371 if (!ifValuesPass && (saveDj * dualIn_ < test1  fabs(saveDj  dualIn_) > checkValue * (1.0 + fabs(saveDj))  fabs(dualIn_) < test2)) { 3372 if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) > 1.0e5)) { 3373 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3374 handler_>message(CLP_PRIMAL_DJ, messages_) 3375 << x << sequenceWithin(sequenceIn_) 3376 << saveDj << dualIn_ 3377 << CoinMessageEol; 3378 if (lastGoodIteration_ != numberIterations_) { 3379 clearAll(); 3380 pivotRow_ = 1; // say no weights update 3381 returnCode = 4; 3382 if (lastGoodIteration_ + 1 == numberIterations_) { 3383 // not looking wonderful  try cleaning bounds 3384 // put nonbasics to bounds in case tolerance moved 3385 abcNonLinearCost_>checkInfeasibilities(0.0); 3386 } 3387 sequenceOut_ = 1; 3388 sequenceIn_ = 1; 3389 if (abcFactorization_>pivots() < 10 && abcFactorization_>pivotTolerance() < 0.25) 3390 abcFactorization_>saferTolerances(1.0, 1.03); 3391 break; 3392 } else { 3393 // take on more relaxed criterion 3394 if (saveDj * dualIn_ < test1  fabs(saveDj  dualIn_) > 2.0e1 * (1.0 + fabs(dualIn_))  fabs(dualIn_) < test2) { 3395 // need to reject something 3396 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3397 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3398 << x << sequenceWithin(sequenceIn_) 3399 << CoinMessageEol; 3400 setFlagged(sequenceIn_); 3401 abcProgress_.incrementTimesFlagged(); 3402 #if 1 //def FEB_TRY \ 3403 // Make safer? 3404 double tolerance = abcFactorization_>pivotTolerance(); 3405 abcFactorization_>saferTolerances(1.0, 1.03); 3406 #endif 3407 abcProgress_.clearBadTimes(); 3408 lastBadIteration_ = numberIterations_; // say be more cautious 3409 clearAll(); 3410 pivotRow_ = 1; 3411 returnCode = 5; 3412 if (tolerance < abcFactorization_>pivotTolerance()) 3413 returnCode = 4; 3414 sequenceOut_ = 1; 3415 sequenceIn_ = 1; 3416 break; 3417 } 3418 } 3461 3419 } else { 3462 3420 //printf("%d %g %g\n",numberIterations_,saveDj,dualIn_); 3463 3421 } 3464 3422 } … … 3472 3430 //abcFactorization_>checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_); 3473 3431 //usefulArray_[arrayForReplaceColumn_].print(); 3474 ftAlpha_ =abcFactorization_>checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],pivotRow_);3475 int updateStatus =abcFactorization_>checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);3432 ftAlpha_ = abcFactorization_>checkReplacePart1(&usefulArray_[arrayForReplaceColumn_], pivotRow_); 3433 int updateStatus = abcFactorization_>checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_); 3476 3434 abcFactorization_>replaceColumnPart3(this, 3477 3478 3479 3480 ftAlpha_?ftAlpha_:alpha_);3481 3435 &usefulArray_[arrayForReplaceColumn_], 3436 &usefulArray_[arrayForFtran_], 3437 pivotRow_, 3438 ftAlpha_ ? ftAlpha_ : alpha_); 3439 3482 3440 // if no pivots, bad update but reasonable alpha  take and invert 3483 if (updateStatus == 2 && 3484 lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e5) 3485 updateStatus = 4; 3441 if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e5) 3442 updateStatus = 4; 3486 3443 if (updateStatus == 1  updateStatus == 4) { 3487 3488 3489 3490 3444 // slight error 3445 if (abcFactorization_>pivots() > 5  updateStatus == 4) { 3446 returnCode = 3; 3447 } 3491 3448 } else if (updateStatus == 2) { 3492 // major error 3493 // better to have small tolerance even if slower 3494 abcFactorization_>zeroTolerance(CoinMin(abcFactorization_>zeroTolerance(), 1.0e15)); 3495 int maxFactor = abcFactorization_>maximumPivots(); 3496 if (maxFactor > 10) { 3497 if (forceFactorization_ < 0) 3498 forceFactorization_ = maxFactor; 3499 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1)); 3500 } 3501 // later we may need to unwind more e.g. fake bounds 3502 if(lastGoodIteration_ != numberIterations_) { 3503 clearAll(); 3504 pivotRow_ = 1; 3505 sequenceIn_ = 1; 3506 sequenceOut_ = 1; 3507 returnCode = 4; 3508 break; 3509 } else { 3510 // need to reject something 3511 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3512 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3513 << x << sequenceWithin(sequenceIn_) 3514 << CoinMessageEol; 3515 setFlagged(sequenceIn_); 3516 abcProgress_.incrementTimesFlagged(); 3517 abcProgress_.clearBadTimes(); 3518 lastBadIteration_ = numberIterations_; // say be more cautious 3519 clearAll(); 3520 pivotRow_ = 1; 3521 sequenceIn_ = 1; 3522 sequenceOut_ = 1; 3523 returnCode = 5; 3524 break; 3525 3526 } 3449 // major error 3450 // better to have small tolerance even if slower 3451 abcFactorization_>zeroTolerance(CoinMin(abcFactorization_>zeroTolerance(), 1.0e15)); 3452 int maxFactor = abcFactorization_>maximumPivots(); 3453 if (maxFactor > 10) { 3454 if (forceFactorization_ < 0) 3455 forceFactorization_ = maxFactor; 3456 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1)); 3457 } 3458 // later we may need to unwind more e.g. fake bounds 3459 if (lastGoodIteration_ != numberIterations_) { 3460 clearAll(); 3461 pivotRow_ = 1; 3462 sequenceIn_ = 1; 3463 sequenceOut_ = 1; 3464 returnCode = 4; 3465 break; 3466 } else { 3467 // need to reject something 3468 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3469 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3470 << x << sequenceWithin(sequenceIn_) 3471 << CoinMessageEol; 3472 setFlagged(sequenceIn_); 3473 abcProgress_.incrementTimesFlagged(); 3474 abcProgress_.clearBadTimes(); 3475 lastBadIteration_ = numberIterations_; // say be more cautious 3476 clearAll(); 3477 pivotRow_ = 1; 3478 sequenceIn_ = 1; 3479 sequenceOut_ = 1; 3480 returnCode = 5; 3481 break; 3482 } 3527 3483 } else if (updateStatus == 3) { 3528 // out of memory 3529 // increase space if not many iterations 3530 if (abcFactorization_>pivots() < 3531 0.5 * abcFactorization_>maximumPivots() && 3532 abcFactorization_>pivots() < 200) 3533 abcFactorization_>areaFactor( 3534 abcFactorization_>areaFactor() * 1.1); 3535 returnCode = 2; // factorize now 3484 // out of memory 3485 // increase space if not many iterations 3486 if (abcFactorization_>pivots() < 0.5 * abcFactorization_>maximumPivots() && abcFactorization_>pivots() < 200) 3487 abcFactorization_>areaFactor( 3488 abcFactorization_>areaFactor() * 1.1); 3489 returnCode = 2; // factorize now 3536 3490 } else if (updateStatus == 5) { 3537 3491 problemStatus_ = 2; // factorize now 3538 3492 } 3539 3493 // here do part of steepest  ready for next iteration 3540 3494 if (!ifValuesPass) 3541 3495 abcPrimalColumnPivot_>updateWeights(&usefulArray_[arrayForFtran_]); 3542 3496 } else { 3543 3497 if (pivotRow_ == 1) { 3544 3545 3546 3547 3548 3549 3550 3498 // no outgoing row is valid 3499 if (valueOut_ != COIN_DBL_MAX) { 3500 double objectiveChange = 0.0; 3501 theta_ = valueOut_  valueIn_; 3502 updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass); 3503 abcSolution_[sequenceIn_] += theta_; 3504 } 3551 3505 #ifdef CLP_USER_DRIVEN1 3552 3506 /* Note if valueOut_ < COIN_DBL_MAX and 3553 3507 theta_ reasonable then this may be a valid sub flip */ 3554 if(!userChoiceValid2(this)) {3555 if (abcFactorization_>pivots()<5) {3556 3557 3558 3559 3560 3561 3562 3563 3564 3565 3566 3567 3568 3569 3570 3571 3572 #endif 3573 3574 3575 3576 3577 3578 3579 3580 3581 3508 if (!userChoiceValid2(this)) { 3509 if (abcFactorization_>pivots() < 5) { 3510 // flag variable 3511 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3512 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3513 << x << sequenceWithin(sequenceIn_) 3514 << CoinMessageEol; 3515 setFlagged(sequenceIn_); 3516 abcProgress_.incrementTimesFlagged(); 3517 abcProgress_.clearBadTimes(); 3518 roundAgain = true; 3519 continue; 3520 } else { 3521 // try refactorizing first 3522 returnCode = 4; //say looks odd but has iterated 3523 break; 3524 } 3525 } 3526 #endif 3527 if (!abcFactorization_>pivots() && acceptablePivot_ <= 1.0e8) { 3528 returnCode = 2; //say looks unbounded 3529 // do ray 3530 primalRay(&usefulArray_[arrayForFtran_]); 3531 } else { 3532 acceptablePivot_ = 1.0e8; 3533 returnCode = 4; //say looks unbounded but has iterated 3534 } 3535 break; 3582 3536 } else { 3583 3584 } 3585 } 3586 3537 // flipping from bound to bound 3538 } 3539 } 3540 3587 3541 double oldCost = 0.0; 3588 3542 if (sequenceOut_ >= 0) 3589 3543 oldCost = abcCost_[sequenceOut_]; 3590 3544 // update primal solution 3591 3545 3592 3546 double objectiveChange = 0.0; 3593 3547 // after this usefulArray_[arrayForFtran_] is not empty  used to update djs 3594 3548 #ifdef CLP_USER_DRIVEN 3595 if (theta_ <0.0) {3596 if (theta_ >=1.0e12)3597 theta_=0.0;3549 if (theta_ < 0.0) { 3550 if (theta_ >= 1.0e12) 3551 theta_ = 0.0; 3598 3552 //else 3599 3553 //printf("negative theta %g\n",theta_); … … 3601 3555 #endif 3602 3556 updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass); 3603 3557 3604 3558 double oldValue = valueIn_; 3605 3559 if (directionIn_ == 1) { 3606 3560 // as if from upper bound 3607 3561 if (sequenceIn_ != sequenceOut_) { 3608 3609 3562 // variable becoming basic 3563 valueIn_ = fabs(theta_); 3610 3564 } else { 3611 3565 valueIn_ = lowerIn_; 3612 3566 } 3613 3567 } else { 3614 3568 // as if from lower bound 3615 3569 if (sequenceIn_ != sequenceOut_) { 3616 3617 3570 // variable becoming basic 3571 valueIn_ += fabs(theta_); 3618 3572 } else { 3619 3573 valueIn_ = upperIn_; 3620 3574 } 3621 3575 } … … 3624 3578 if (sequenceIn_ != sequenceOut_) { 3625 3579 if (directionOut_ > 0) { 3626 3580 valueOut_ = lowerOut_; 3627 3581 } else { 3628 3629 } 3630 if (valueOut_ < abcLower_[sequenceOut_]  primalTolerance_)3631 3582 valueOut_ = upperOut_; 3583 } 3584 if (valueOut_ < abcLower_[sequenceOut_]  primalTolerance_) 3585 valueOut_ = abcLower_[sequenceOut_]  0.9 * primalTolerance_; 3632 3586 else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_) 3633 3587 valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_; 3634 3588 // may not be exactly at bound and bounds may have changed 3635 3589 // Make sure outgoing looks feasible … … 3651 3605 // maximum iterations or equivalent 3652 3606 returnCode = 3; 3653 } else if(numberIterations_ == lastGoodIteration_ 3654 + 2 * abcFactorization_>maximumPivots()) { 3607 } else if (numberIterations_ == lastGoodIteration_ + 2 * abcFactorization_>maximumPivots()) { 3655 3608 // done a lot of flips  be safe 3656 3609 returnCode = 2; // refactorize … … 3660 3613 int status = eventHandler_>event(ClpEventHandler::endOfIteration); 3661 3614 if (status >= 0) { 3662 3663 3664 3615 problemStatus_ = 5; 3616 secondaryStatus_ = ClpEventHandler::endOfIteration; 3617 returnCode = 3; 3665 3618 } 3666 3619 } … … 3678 3631 +3 max iterations (iteration done) 3679 3632 */ 3680 int 3681 AbcSimplexPrimal::pivotResult4(int ifValuesPass) 3633 int AbcSimplexPrimal::pivotResult4(int ifValuesPass) 3682 3634 { 3683 3635 3684 3636 int returnCode = 1; 3685 int numberMinor =0;3686 int numberDone =0;3687 CoinIndexedVector * 3637 int numberMinor = 0; 3638 int numberDone = 0; 3639 CoinIndexedVector *vector[16]; 3688 3640 pivotStruct stuff[4]; 3689 vector[0] =&usefulArray_[arrayForFtran_];3690 vector[1] =&usefulArray_[arrayForFlipBounds_];3691 vector[2] =&usefulArray_[arrayForTableauRow_];3692 vector[3] =&usefulArray_[arrayForBtran_];3641 vector[0] = &usefulArray_[arrayForFtran_]; 3642 vector[1] = &usefulArray_[arrayForFlipBounds_]; 3643 vector[2] = &usefulArray_[arrayForTableauRow_]; 3644 vector[3] = &usefulArray_[arrayForBtran_]; 3693 3645 /* 3694 3646 For pricing we need to get a vector with difference in costs … … 3699 3651 //double * saveCosts=CoinCopyOfArray(costBasic_,numberRows_); 3700 3652 //double saveCostsIn[4]; 3701 for (int iMinor =0;iMinor<4;iMinor++) {3702 int sequenceIn =multipleSequenceIn_[iMinor];3703 if (sequenceIn <0)3653 for (int iMinor = 0; iMinor < 4; iMinor++) { 3654 int sequenceIn = multipleSequenceIn_[iMinor]; 3655 if (sequenceIn < 0) 3704 3656 break; 3705 stuff[iMinor].valuesPass_ =ifValuesPass;3706 stuff[iMinor].lowerIn_ =abcLower_[sequenceIn];3707 stuff[iMinor].upperIn_ =abcUpper_[sequenceIn];3708 stuff[iMinor].valueIn_ =abcSolution_[sequenceIn];3709 stuff[iMinor].sequenceIn_ =sequenceIn;3657 stuff[iMinor].valuesPass_ = ifValuesPass; 3658 stuff[iMinor].lowerIn_ = abcLower_[sequenceIn]; 3659 stuff[iMinor].upperIn_ = abcUpper_[sequenceIn]; 3660 stuff[iMinor].valueIn_ = abcSolution_[sequenceIn]; 3661 stuff[iMinor].sequenceIn_ = sequenceIn; 3710 3662 //saveCostsIn[iMinor]=abcCost_[sequenceIn]; 3711 3663 numberMinor++; 3712 3664 if (iMinor) { 3713 vector[4 *iMinor]=rowArray_[2*iMinor2];3714 vector[4 *iMinor+1]=rowArray_[2*iMinor1];3715 vector[4 *iMinor+2]=columnArray_[2*iMinor2];3716 vector[4 *iMinor+3]=columnArray_[2*iMinor1];3717 } 3718 for (int i =0;i<4;i++)3719 vector[4 *iMinor+i]>checkClear();3720 unpack(*vector[4 *iMinor],sequenceIn);3721 } 3722 int numberLeft =numberMinor;3665 vector[4 * iMinor] = rowArray_[2 * iMinor  2]; 3666 vector[4 * iMinor + 1] = rowArray_[2 * iMinor  1]; 3667 vector[4 * iMinor + 2] = columnArray_[2 * iMinor  2]; 3668 vector[4 * iMinor + 3] = columnArray_[2 * iMinor  1]; 3669 } 3670 for (int i = 0; i < 4; i++) 3671 vector[4 * iMinor + i]>checkClear(); 3672 unpack(*vector[4 * iMinor], sequenceIn); 3673 } 3674 int numberLeft = numberMinor; 3723 3675 // parallel (with cpu) 3724 for (int iMinor =1;iMinor<numberMinor;iMinor++) {3676 for (int iMinor = 1; iMinor < numberMinor; iMinor++) { 3725 3677 // update the incoming columns 3726 cilk_spawn abcFactorization_>updateColumnFT(*vector[4 *iMinor],*vector[4*iMinor+3],iMinor);3727 } 3728 abcFactorization_>updateColumnFT(*vector[0], *vector[+3],0);3678 cilk_spawn abcFactorization_>updateColumnFT(*vector[4 * iMinor], *vector[4 * iMinor + 3], iMinor); 3679 } 3680 abcFactorization_>updateColumnFT(*vector[0], *vector[+3], 0); 3729 3681 cilk_sync; 3730 for (int iMinor =0;iMinor<numberMinor;iMinor++) {3682 for (int iMinor = 0; iMinor < numberMinor; iMinor++) { 3731 3683 // find best (or first if values pass) 3732 int numberDo =1;3733 int jMinor =0;3734 while (jMinor <numberDo) {3735 int sequenceIn =stuff[jMinor].sequenceIn_;3736 double dj =abcDj_[sequenceIn];3737 bool bad =false;3684 int numberDo = 1; 3685 int jMinor = 0; 3686 while (jMinor < numberDo) { 3687 int sequenceIn = stuff[jMinor].sequenceIn_; 3688 double dj = abcDj_[sequenceIn]; 3689 bool bad = false; 3738 3690 if (!bad) { 3739 stuff[jMinor].dualIn_=dj;3740 stuff[jMinor].saveDualIn_=dj;3741 if (dj<0.0)3742 stuff[jMinor].directionIn_=1;3743 3744 stuff[jMinor].directionIn_=1;3745 3691 stuff[jMinor].dualIn_ = dj; 3692 stuff[jMinor].saveDualIn_ = dj; 3693 if (dj < 0.0) 3694 stuff[jMinor].directionIn_ = 1; 3695 else 3696 stuff[jMinor].directionIn_ = 1; 3697 jMinor++; 3746 3698 } else { 3747 3748 3749 3750 for (int i=0;i<4;i++) {3751 vector[4*jMinor+i]>clear();3752 vector[4*jMinor+i]=vector[4*numberDo+i];3753 vector[4*numberDo+i]=NULL;3754 3755 stuff[jMinor]=stuff[numberDo];3756 } 3757 } 3758 for (int jMinor =0;jMinor<numberDo;jMinor++) {3699 numberDo; 3700 numberLeft; 3701 // throw away 3702 for (int i = 0; i < 4; i++) { 3703 vector[4 * jMinor + i]>clear(); 3704 vector[4 * jMinor + i] = vector[4 * numberDo + i]; 3705 vector[4 * numberDo + i] = NULL; 3706 } 3707 stuff[jMinor] = stuff[numberDo]; 3708 } 3709 } 3710 for (int jMinor = 0; jMinor < numberDo; jMinor++) { 3759 3711 // do ratio test and recompute dj 3760 primalRow(vector[4 *jMinor],vector[4*jMinor+1],vector[4*jMinor+2],stuff[jMinor]);3712 primalRow(vector[4 * jMinor], vector[4 * jMinor + 1], vector[4 * jMinor + 2], stuff[jMinor]); 3761 3713 } 3762 3714 // choose best 3763 int iBest =1;3764 double bestMovement =COIN_DBL_MAX;3765 for (int jMinor =0;jMinor<numberDo;jMinor++) {3766 double movement =stuff[jMinor].theta_*fabs(stuff[jMinor].dualIn_);3767 if (movement >bestMovement) {3768 bestMovement=movement;3769 iBest=jMinor;3715 int iBest = 1; 3716 double bestMovement = COIN_DBL_MAX; 3717 for (int jMinor = 0; jMinor < numberDo; jMinor++) { 3718 double movement = stuff[jMinor].theta_ * fabs(stuff[jMinor].dualIn_); 3719 if (movement > bestMovement) { 3720 bestMovement = movement; 3721 iBest = jMinor; 3770 3722 } 3771 3723 } … … 3774 3726 iBest=0; 3775 3727 #endif 3776 if (iBest >=0) {3777 dualIn_ =stuff[iBest].dualIn_;3778 dualOut_ =stuff[iBest].dualOut_;3779 lowerOut_ =stuff[iBest].lowerOut_;3780 upperOut_ =stuff[iBest].upperOut_;3781 valueOut_ =stuff[iBest].valueOut_;3782 sequenceOut_ =stuff[iBest].sequenceOut_;3783 sequenceIn_ =stuff[iBest].sequenceIn_;3784 lowerIn_ =stuff[iBest].lowerIn_;3785 upperIn_ =stuff[iBest].upperIn_;3786 valueIn_ =stuff[iBest].valueIn_;3787 directionIn_ =stuff[iBest].directionIn_;3788 directionOut_ =stuff[iBest].directionOut_;3789 pivotRow_ =stuff[iBest].pivotRow_;3790 theta_ =stuff[iBest].theta_;3791 alpha_ =stuff[iBest].alpha_;3728 if (iBest >= 0) { 3729 dualIn_ = stuff[iBest].dualIn_; 3730 dualOut_ = stuff[iBest].dualOut_; 3731 lowerOut_ = stuff[iBest].lowerOut_; 3732 upperOut_ = stuff[iBest].upperOut_; 3733 valueOut_ = stuff[iBest].valueOut_; 3734 sequenceOut_ = stuff[iBest].sequenceOut_; 3735 sequenceIn_ = stuff[iBest].sequenceIn_; 3736 lowerIn_ = stuff[iBest].lowerIn_; 3737 upperIn_ = stuff[iBest].upperIn_; 3738 valueIn_ = stuff[iBest].valueIn_; 3739 directionIn_ = stuff[iBest].directionIn_; 3740 directionOut_ = stuff[iBest].directionOut_; 3741 pivotRow_ = stuff[iBest].pivotRow_; 3742 theta_ = stuff[iBest].theta_; 3743 alpha_ = stuff[iBest].alpha_; 3792 3744 #ifdef MULTIPLE_PRICE 3793 for (int i =0;i<4*numberLeft;i++)3794 3745 for (int i = 0; i < 4 * numberLeft; i++) 3746 vector[i]>clear(); 3795 3747 return 0; 3796 3748 #endif … … 3799 3751 double theta2 = numberIterations_ * abcNonLinearCost_>averageTheta(); 3800 3752 // Set average theta 3801 abcNonLinearCost_>setAverageTheta((theta1 + theta2) / (static_cast< double>(numberIterations_ + 1)));3753 abcNonLinearCost_>setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1))); 3802 3754 if (pivotRow_ == 1  (pivotRow_ >= 0 && fabs(alpha_) < 1.0e5)) { 3803 if(fabs(dualIn_) < 1.0e2 * dualTolerance_) {3804 3805 3806 stuff[iBest].valuesPass_=0;3807 primalRow(vector[4*iBest],vector[4*iBest+1],vector[4*iBest+2],stuff[iBest]);3808 dualIn_=stuff[iBest].dualIn_;3809 dualOut_=stuff[iBest].dualOut_;3810 lowerOut_=stuff[iBest].lowerOut_;3811 upperOut_=stuff[iBest].upperOut_;3812 valueOut_=stuff[iBest].valueOut_;3813 sequenceOut_=stuff[iBest].sequenceOut_;3814 sequenceIn_=stuff[iBest].sequenceIn_;3815 lowerIn_=stuff[iBest].lowerIn_;3816 upperIn_=stuff[iBest].upperIn_;3817 valueIn_=stuff[iBest].valueIn_;3818 directionIn_=stuff[iBest].directionIn_;3819 directionOut_=stuff[iBest].directionOut_;3820 pivotRow_=stuff[iBest].pivotRow_;3821 theta_=stuff[iBest].theta_;3822 alpha_=stuff[iBest].alpha_;3823 3824 3825 3826 3827 3828 3829 3830 3831 3832 3833 3834 3835 3836 3837 3838 3839 } 3840 CoinIndexedVector * bestUpdate=NULL;3755 if (fabs(dualIn_) < 1.0e2 * dualTolerance_) { 3756 // try other way 3757 stuff[iBest].directionIn_ = directionIn_; 3758 stuff[iBest].valuesPass_ = 0; 3759 primalRow(vector[4 * iBest], vector[4 * iBest + 1], vector[4 * iBest + 2], stuff[iBest]); 3760 dualIn_ = stuff[iBest].dualIn_; 3761 dualOut_ = stuff[iBest].dualOut_; 3762 lowerOut_ = stuff[iBest].lowerOut_; 3763 upperOut_ = stuff[iBest].upperOut_; 3764 valueOut_ = stuff[iBest].valueOut_; 3765 sequenceOut_ = stuff[iBest].sequenceOut_; 3766 sequenceIn_ = stuff[iBest].sequenceIn_; 3767 lowerIn_ = stuff[iBest].lowerIn_; 3768 upperIn_ = stuff[iBest].upperIn_; 3769 valueIn_ = stuff[iBest].valueIn_; 3770 directionIn_ = stuff[iBest].directionIn_; 3771 directionOut_ = stuff[iBest].directionOut_; 3772 pivotRow_ = stuff[iBest].pivotRow_; 3773 theta_ = stuff[iBest].theta_; 3774 alpha_ = stuff[iBest].alpha_; 3775 } 3776 if (pivotRow_ == 1  (pivotRow_ >= 0 && fabs(alpha_) < 1.0e5)) { 3777 // reject it 3778 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 3779 handler_>message(CLP_SIMPLEX_FLAG, messages_) 3780 << x << sequenceWithin(sequenceIn_) 3781 << CoinMessageEol; 3782 setFlagged(sequenceIn_); 3783 abcProgress_.incrementTimesFlagged(); 3784 abcProgress_.clearBadTimes(); 3785 lastBadIteration_ = numberIterations_; // say be more cautious 3786 clearAll(); 3787 pivotRow_ = 1; 3788 returnCode = 5; 3789 break; 3790 } 3791 } 3792 CoinIndexedVector *bestUpdate = NULL; 3841 3793 if (pivotRow_ >= 0) { 3842 3843 3844 3845 CoinIndexedVector *tempVector[4];3846 for (int i=0;i<4;i++) 3847 tempVector[i] = vector[4*iBest+i];3848 3849 bestUpdate=tempVector[0];3850 3851 3852 3853 3854 3855 for (int i=0;i<4;i++) {3856 vector[4*iBest+i]=vector[4*numberLeft+i];3857 vector[4*numberLeft+i]=NULL;3858 3859 stuff[iBest]=stuff[numberLeft];3860 3861 3862 for (int jMinor=0;jMinor<numberLeft;jMinor++) {3863 cilk_spawn updatePartialUpdate(*vector[4*jMinor+3]);3864 3865 3866 for (int jMinor=0;jMinor<numberLeft;jMinor++) {3867 stuff[jMinor].dualIn_=cilk_spawn updateMinorCandidate(*bestUpdate,*vector[4*jMinor],stuff[jMinor].sequenceIn_);3868 3869 3870 3871 for (int i=1;i<4;i++) 3872 3873 if (returnCode<3) {3874 3875 3876 3877 3794 // Normal pivot 3795 // if stable replace in basis 3796 // check update 3797 CoinIndexedVector *tempVector[4]; 3798 for (int i = 0; i < 4; i++) 3799 tempVector[i] = vector[4 * iBest + i]; 3800 returnCode = cilk_spawn doFTUpdate(tempVector); 3801 bestUpdate = tempVector[0]; 3802 // after this bestUpdate is not empty  used to update djs ?? 3803 cilk_spawn updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass); 3804 numberDone++; 3805 numberLeft; 3806 // throw away 3807 for (int i = 0; i < 4; i++) { 3808 vector[4 * iBest + i] = vector[4 * numberLeft + i]; 3809 vector[4 * numberLeft + i] = NULL; 3810 } 3811 stuff[iBest] = stuff[numberLeft]; 3812 // update pi and other vectors and FT stuff 3813 // parallel (can go 8 way?) 3814 for (int jMinor = 0; jMinor < numberLeft; jMinor++) { 3815 cilk_spawn updatePartialUpdate(*vector[4 * jMinor + 3]); 3816 } 3817 // parallel 3818 for (int jMinor = 0; jMinor < numberLeft; jMinor++) { 3819 stuff[jMinor].dualIn_ = cilk_spawn updateMinorCandidate(*bestUpdate, *vector[4 * jMinor], stuff[jMinor].sequenceIn_); 3820 } 3821 cilk_sync; 3822 // throw away 3823 for (int i = 1; i < 4; i++) 3824 tempVector[i]>clear(); 3825 if (returnCode < 3) { 3826 clearAll(); 3827 break; 3828 } 3829 // end Normal pivot 3878 3830 } else { 3879 3880 vector[4*iBest+3]>clear();3881 3882 3883 3884 3885 updatePrimalsInPrimal(*vector[4*iBest], theta_, ifValuesPass);3886 3887 3888 3889 3890 3891 primalRay(vector[4*iBest]);3892 3893 3894 3895 3896 3897 3898 3899 bestUpdate=vector[4*iBest];3900 3901 3902 3903 3904 3905 3906 for (int i=0;i<4;i++) {3907 3908 vector[4*iBest+i]>clear();3909 vector[4*iBest+i]=vector[4*numberLeft+i];3910 vector[4*numberLeft+i]=NULL;3911 3912 stuff[iBest]=stuff[numberLeft];3913 3914 for (int jMinor=0;jMinor<numberLeft;jMinor++) {3915 stuff[jMinor].dualIn_=cilk_spawn updateMinorCandidate(*bestUpdate,*vector[4*jMinor],stuff[jMinor].sequenceIn_);3916 3917 3918 3919 3831 // start Flip 3832 vector[4 * iBest + 3]>clear(); 3833 if (pivotRow_ == 1) { 3834 // no outgoing row is valid 3835 if (valueOut_ != COIN_DBL_MAX) { 3836 theta_ = valueOut_  valueIn_; 3837 updatePrimalsInPrimal(*vector[4 * iBest], theta_, ifValuesPass); 3838 abcSolution_[sequenceIn_] += theta_; 3839 } 3840 if (!abcFactorization_>pivots() && acceptablePivot_ <= 1.0e8) { 3841 returnCode = 2; //say looks unbounded 3842 // do ray 3843 primalRay(vector[4 * iBest]); 3844 } else { 3845 acceptablePivot_ = 1.0e8; 3846 returnCode = 4; //say looks unbounded but has iterated 3847 } 3848 break; 3849 } else { 3850 // flipping from bound to bound 3851 bestUpdate = vector[4 * iBest]; 3852 // after this bestUpdate is not empty  used to update djs ?? 3853 updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass); 3854 // throw away best BUT remember we are updating costs somehow 3855 numberDone++; 3856 numberLeft; 3857 // throw away 3858 for (int i = 0; i < 4; i++) { 3859 if (i) 3860 vector[4 * iBest + i]>clear(); 3861 vector[4 * iBest + i] = vector[4 * numberLeft + i]; 3862 vector[4 * numberLeft + i] = NULL; 3863 } 3864 stuff[iBest] = stuff[numberLeft]; 3865 // parallel 3866 for (int jMinor = 0; jMinor < numberLeft; jMinor++) { 3867 stuff[jMinor].dualIn_ = cilk_spawn updateMinorCandidate(*bestUpdate, *vector[4 * jMinor], stuff[jMinor].sequenceIn_); 3868 } 3869 cilk_sync; 3870 } 3871 // end Flip 3920 3872 } 3921 3873 bestUpdate>clear(); // as only works in values pass 3922 3874 3923 3875 if (directionIn_ == 1) { 3924 3925 3926 3927 3928 3929 3930 3876 // as if from upper bound 3877 if (sequenceIn_ != sequenceOut_) { 3878 // variable becoming basic 3879 valueIn_ = fabs(theta_); 3880 } else { 3881 valueIn_ = lowerIn_; 3882 } 3931 3883 } else { 3932 3933 3934 3935 3936 3937 3938 3884 // as if from lower bound 3885 if (sequenceIn_ != sequenceOut_) { 3886 // variable becoming basic 3887 valueIn_ += fabs(theta_); 3888 } else { 3889 valueIn_ = upperIn_; 3890 } 3939 3891 } 3940 3892 // outgoing 3941 3893 if (sequenceIn_ != sequenceOut_) { 3942 3943 3944 3945 3946 3947 if(valueOut_ < abcLower_[sequenceOut_]  primalTolerance_)3948 3949 3950 3951 3952 3953 3954 3894 if (directionOut_ > 0) { 3895 valueOut_ = lowerOut_; 3896 } else { 3897 valueOut_ = upperOut_; 3898 } 3899 if (valueOut_ < abcLower_[sequenceOut_]  primalTolerance_) 3900 valueOut_ = abcLower_[sequenceOut_]  0.9 * primalTolerance_; 3901 else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_) 3902 valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_; 3903 // may not be exactly at bound and bounds may have changed 3904 // Make sure outgoing looks feasible 3905 directionOut_ = abcNonLinearCost_>setOneOutgoing(pivotRow_, valueOut_); 3906 abcSolution_[sequenceOut_] = valueOut_; 3955 3907 } 3956 3908 // change cost and bounds on incoming if primal … … 3959 3911 swapPrimalStuff(); 3960 3912 if (whatNext == 1) { 3961 3913 returnCode = 2; // refactorize 3962 3914 } else if (whatNext == 2) { 3963 // maximum iterations or equivalent 3964 returnCode = 3; 3965 } else if(numberIterations_ == lastGoodIteration_ 3966 + 2 * abcFactorization_>maximumPivots()) { 3967 // done a lot of flips  be safe 3968 returnCode = 2; // refactorize 3915 // maximum iterations or equivalent 3916 returnCode = 3; 3917 } else if (numberIterations_ == lastGoodIteration_ + 2 * abcFactorization_>maximumPivots()) { 3918 // done a lot of flips  be safe 3919 returnCode = 2; // refactorize 3969 3920 } 3970 3921 // Check event 3971 3922 { 3972 3973 3974 3975 3976 3977 3978 } 3979 } 3980 if (returnCode !=1)3923 int status = eventHandler_>event(ClpEventHandler::endOfIteration); 3924 if (status >= 0) { 3925 problemStatus_ = 5; 3926 secondaryStatus_ = ClpEventHandler::endOfIteration; 3927 returnCode = 3; 3928 } 3929 } 3930 } 3931 if (returnCode != 1) 3981 3932 break; 3982 3933 } 3983 for (int i =0;i<16;i++) {3934 for (int i = 0; i < 16; i++) { 3984 3935 if (vector[i]) 3985 3936 vector[i]>clear(); … … 3991 3942 } 3992 3943 // Do FT update as separate function for minor iterations (nonzero return code on problems) 3993 int 3994 AbcSimplexPrimal::doFTUpdate(CoinIndexedVector * vector[4]) 3944 int AbcSimplexPrimal::doFTUpdate(CoinIndexedVector *vector[4]) 3995 3945 { 3996 ftAlpha_ =abcFactorization_>checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],3997 vector[3],pivotRow_);3998 int updateStatus =abcFactorization_>checkReplacePart2(pivotRow_,btranAlpha_,alpha_,ftAlpha_);3946 ftAlpha_ = abcFactorization_>checkReplacePart1(&usefulArray_[arrayForReplaceColumn_], 3947 vector[3], pivotRow_); 3948 int updateStatus = abcFactorization_>checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_); 3999 3949 if (!updateStatus) 4000 3950 abcFactorization_>replaceColumnPart3(this, 4001 4002 4003 4004 4005 ftAlpha_?ftAlpha_:alpha_);3951 &usefulArray_[arrayForReplaceColumn_], 3952 vector[0], 3953 vector[3], 3954 pivotRow_, 3955 ftAlpha_ ? ftAlpha_ : alpha_); 4006 3956 else 4007 3957 vector[3]>clear(); 4008 int returnCode =0;3958 int returnCode = 0; 4009 3959 // if no pivots, bad update but reasonable alpha  take and invert 4010 if (updateStatus == 2 && 4011 lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e5) 3960 if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e5) 4012 3961 updateStatus = 4; 4013 3962 if (updateStatus == 1  updateStatus == 4) { … … 4023 3972 if (maxFactor > 10) { 4024 3973 if (forceFactorization_ < 0) 4025 3974 forceFactorization_ = maxFactor; 4026 3975 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1)); 4027 3976 } 4028 3977 // later we may need to unwind more e.g. fake bounds 4029 if (lastGoodIteration_ != numberIterations_) {3978 if (lastGoodIteration_ != numberIterations_) { 4030 3979 pivotRow_ = 1; 4031 3980 sequenceIn_ = 1; … … 4036 3985 char x = isColumn(sequenceIn_) ? 'C' : 'R'; 4037 3986 handler_>message(CLP_SIMPLEX_FLAG, messages_) 4038 4039 3987 << x << sequenceWithin(sequenceIn_) 3988 << CoinMessageEol; 4040 3989 setFlagged(sequenceIn_); 4041 3990 abcProgress_.incrementTimesFlagged(); … … 4046 3995 sequenceOut_ = 1; 4047 3996 returnCode = 5; 4048 4049 3997 } 4050 3998 } else if (updateStatus == 3) { 4051 3999 // out of memory 4052 4000 // increase space if not many iterations 4053 if (abcFactorization_>pivots() < 4054 0.5 * abcFactorization_>maximumPivots() && 4055 abcFactorization_>pivots() < 200) 4001 if (abcFactorization_>pivots() < 0.5 * abcFactorization_>maximumPivots() && abcFactorization_>pivots() < 200) 4056 4002 abcFactorization_>areaFactor( 4057 4003 abcFactorization_>areaFactor() * 1.1); 4058 4004 returnCode = 2; // factorize now 4059 4005 } else if (updateStatus == 5) { 4060 4006 problemStatus_ = 2; // factorize now 4061 returnCode =2;4007 returnCode = 2; 4062 4008 } 4063 4009 return returnCode; … … 4066 4012 costs are changed 4067 4013 */ 4068 void 4069 AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector & rowArray, 4070 double theta,bool valuesPass) 4014 void AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector &rowArray, 4015 double theta, bool valuesPass) 4071 4016 { 4072 4017 // Cost on pivot row may change  may need to change dualIn … … 4074 4019 //if (pivotRow_ >= 0) 4075 4020 //oldCost = abcCost_[sequenceOut_]; 4076 double * 4021 double *work = rowArray.denseVector(); 4077 4022 int number = rowArray.getNumElements(); 4078 int * 4023 int *which = rowArray.getIndices(); 4079 4024 // allow for case where bound+tolerance == bound 4080 4025 //double tolerance = 0.999999*primalTolerance_; … … 4083 4028 if (!valuesPass) { 4084 4029 for (iIndex = 0; iIndex < number; iIndex++) { 4085 4030 4086 4031 int iRow = which[iIndex]; 4087 4032 double alpha = work[iRow]; … … 4094 4039 solutionBasic_[iRow] = value; 4095 4040 if (active(iRow)  theta_ < 0.0) { 4096 4097 4098 4099 4100 4101 4102 4103 4104 4105 4106 4107 4108 4109 4110 4111 4112 4041 clearActive(iRow); 4042 // But make sure one going out is feasible 4043 if (change > 0.0) { 4044 // going down 4045 if (value <= abcLower_[iPivot] + primalTolerance_) { 4046 if (iPivot == sequenceOut_ && value >= abcLower_[iPivot]  relaxedTolerance) 4047 value = abcLower_[iPivot]; 4048 abcNonLinearCost_>setOneBasic(iRow, value); 4049 } 4050 } else { 4051 // going up 4052 if (value >= abcUpper_[iPivot]  primalTolerance_) { 4053 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance) 4054 value = abcUpper_[iPivot]; 4055 abcNonLinearCost_>setOneBasic(iRow, value); 4056 } 4057 } 4113 4058 } 4114 4059 } … … 4116 4061 // values pass so look at all 4117 4062 for (iIndex = 0; iIndex < number; iIndex++) { 4118 4063 4119 4064 int iRow = which[iIndex]; 4120 4065 double alpha = work[iRow]; …