Changeset 212 for branches/pre/ClpPrimalColumnSteepest.cpp
 Timestamp:
 Oct 2, 2003 1:21:02 PM (17 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/pre/ClpPrimalColumnSteepest.cpp
r210 r212 160 160 { 161 161 assert(model_); 162 if (model_>nonLinearCost()>lookBothWays()model_>algorithm()==2) { 163 // Do old way 164 return pivotColumnOldMethod(updates,spareRow1,spareRow2, 165 spareColumn1,spareColumn2); 166 } 167 int number=0; 168 int * index; 169 // dj could be very small (or even zero  take care) 170 double dj = model_>dualIn(); 171 double tolerance=model_>currentDualTolerance(); 172 // we can't really trust infeasibilities if there is dual error 173 // this coding has to mimic coding in checkDualSolution 174 double error = min(1.0e3,model_>largestDualError()); 175 // allow tolerance at least slightly bigger than standard 176 tolerance = tolerance + error; 177 int pivotRow = model_>pivotRow(); 178 int anyUpdates; 179 double * infeas = infeasible_>denseVector(); 180 181 // Local copy of mode so can decide what to do 182 int switchType; 183 if (mode_==4) 184 switchType = 5numberSwitched_; 185 else 186 switchType = mode_; 187 /* switchType  188 0  all exact devex 189 1  all steepest 190 2  some exact devex 191 3  auto some exact devex 192 4  devex 193 5  dantzig 194 */ 195 if (updates>getNumElements()) { 196 // would have to have two goes for devex, three for steepest 197 anyUpdates=2; 198 // add in pivot contribution 199 if (pivotRow>=0) 200 updates>add(pivotRow,dj); 201 } else if (pivotRow>=0) { 202 if (fabs(dj)>1.0e15) { 203 // some dj 204 updates>insert(pivotRow,dj); 205 if (fabs(dj)>1.0e6) { 206 // reasonable size 207 anyUpdates=1; 208 } else { 209 // too small 210 anyUpdates=2; 211 } 212 } else { 213 // zero dj 214 anyUpdates=1; 215 } 216 } else if (pivotSequence_>=0){ 217 // just after refactorization 218 anyUpdates=1; 219 } else { 220 // sub flip  nothing to do 221 anyUpdates=0; 222 } 223 int sequenceOut = model_>sequenceOut(); 224 if (switchType==5) { 225 if (anyUpdates>0) { 226 justDjs(updates,spareRow1,spareRow2, 227 spareColumn1,spareColumn2); 228 } 229 } else if (anyUpdates==1) { 230 if (switchType<4) { 231 // exact etc when can use dj 232 djsAndSteepest(updates,spareRow1,spareRow2, 233 spareColumn1,spareColumn2); 234 } else { 235 // devex etc when can use dj 236 djsAndDevex(updates,spareRow1,spareRow2, 237 spareColumn1,spareColumn2); 238 } 239 } else if (anyUpdates==1) { 240 if (switchType<4) { 241 // exact etc when djs okay 242 justSteepest(updates,spareRow1,spareRow2, 243 spareColumn1,spareColumn2); 244 } else { 245 // devex etc when djs okay 246 justDevex(updates,spareRow1,spareRow2, 247 spareColumn1,spareColumn2); 248 } 249 } else if (anyUpdates==2) { 250 if (switchType<4) { 251 // exact etc when have to use pivot 252 djsAndSteepest2(updates,spareRow1,spareRow2, 253 spareColumn1,spareColumn2); 254 } else { 255 // devex etc when have to use pivot 256 djsAndDevex2(updates,spareRow1,spareRow2, 257 spareColumn1,spareColumn2); 258 } 259 } 260 // make sure outgoing from last iteration okay 261 if (sequenceOut>=0) { 262 ClpSimplex::Status status = model_>getStatus(sequenceOut); 263 double value = model_>reducedCost(sequenceOut); 264 265 switch(status) { 266 267 case ClpSimplex::basic: 268 case ClpSimplex::isFixed: 269 break; 270 case ClpSimplex::isFree: 271 case ClpSimplex::superBasic: 272 if (fabs(value)>FREE_ACCEPT*tolerance) { 273 // we are going to bias towards free (but only if reasonable) 274 value *= FREE_BIAS; 275 // store square in list 276 if (infeas[sequenceOut]) 277 infeas[sequenceOut] = value*value; // already there 278 else 279 infeasible_>quickAdd(sequenceOut,value*value); 280 } else { 281 infeasible_>zero(sequenceOut); 282 } 283 break; 284 case ClpSimplex::atUpperBound: 285 if (value>tolerance) { 286 // store square in list 287 if (infeas[sequenceOut]) 288 infeas[sequenceOut] = value*value; // already there 289 else 290 infeasible_>quickAdd(sequenceOut,value*value); 291 } else { 292 infeasible_>zero(sequenceOut); 293 } 294 break; 295 case ClpSimplex::atLowerBound: 296 if (value<tolerance) { 297 // store square in list 298 if (infeas[sequenceOut]) 299 infeas[sequenceOut] = value*value; // already there 300 else 301 infeasible_>quickAdd(sequenceOut,value*value); 302 } else { 303 infeasible_>zero(sequenceOut); 304 } 305 } 306 } 307 // update of duals finished  now do pricing 308 // See what sort of pricing 309 int numberWanted=0; 310 number = infeasible_>getNumElements(); 311 int numberColumns = model_>numberColumns(); 312 if (switchType==5) { 313 pivotSequence_=1; 314 pivotRow=1; 315 // See if to switch 316 int numberElements = model_>factorization()>numberElements(); 317 int numberRows = model_>numberRows(); 318 // ratio is done on number of columns here 319 //double ratio = (double) numberElements/(double) numberColumns; 320 double ratio = (double) numberElements/(double) numberRows; 321 //double ratio = (double) numberElements/(double) model_>clpMatrix()>getNumElements(); 322 if (ratio<0.1) { 323 numberWanted = max(100,number/200); 324 } else if (ratio<0.15) { 325 numberWanted = max(500,number/40); 326 } else if (ratio<0.2) { 327 numberWanted = max(2000,number/10); 328 numberWanted = max(numberWanted,numberColumns/30); 329 } else { 330 switchType=4; 331 // initialize 332 numberSwitched_++; 333 // Make sure will redo 334 delete [] weights_; 335 weights_=NULL; 336 saveWeights(model_,4); 337 printf("switching to devex %d nel ratio %g\n",numberElements,ratio); 338 } 339 if (model_>numberIterations()%1000==0) 340 printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted); 341 } 342 if(switchType==4) { 343 // Still in devex mode 344 int numberElements = model_>factorization()>numberElements(); 345 int numberRows = model_>numberRows(); 346 // ratio is done on number of rows here 347 double ratio = (double) numberElements/(double) numberRows; 348 // Go to steepest if lot of iterations? 349 if (ratio<1.0) { 350 numberWanted = max(2000,number/20); 351 } else if (ratio<2.0) { 352 numberWanted = max(2000,number/10); 353 numberWanted = max(numberWanted,numberColumns/20); 354 } else { 355 // we can zero out 356 updates>clear(); 357 spareColumn1>clear(); 358 switchType=3; 359 // initialize 360 pivotSequence_=1; 361 pivotRow=1; 362 numberSwitched_++; 363 // Make sure will redo 364 delete [] weights_; 365 weights_=NULL; 366 saveWeights(model_,4); 367 printf("switching to exact %d nel ratio %g\n",numberElements,ratio); 368 updates>clear(); 369 } 370 if (model_>numberIterations()%1000==0) 371 printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted); 372 } 373 if (switchType<4) { 374 if (switchType<2 ) { 375 numberWanted = number+1; 376 } else if (switchType==2) { 377 numberWanted = max(2000,number/8); 378 } else { 379 int numberElements = model_>factorization()>numberElements(); 380 double ratio = (double) numberElements/(double) model_>numberRows(); 381 if (ratio<1.0) { 382 numberWanted = max(2000,number/20); 383 } else if (ratio<5.0) { 384 numberWanted = max(2000,number/10); 385 numberWanted = max(numberWanted,numberColumns/20); 386 } else if (ratio<10.0) { 387 numberWanted = max(2000,number/8); 388 numberWanted = max(numberWanted,numberColumns/20); 389 } else { 390 ratio = number * (ratio/80.0); 391 if (ratio>number) { 392 numberWanted=number+1; 393 } else { 394 numberWanted = max(2000,(int) ratio); 395 numberWanted = max(numberWanted,numberColumns/10); 396 } 397 } 398 } 399 } 400 401 402 double bestDj = 1.0e30; 403 int bestSequence=1; 404 405 int i,iSequence; 406 index = infeasible_>getIndices(); 407 number = infeasible_>getNumElements(); 408 // Resort infeasible every 100 pivots 409 if (0&&model_>factorization()>pivots()>0&& 410 (model_>factorization()>pivots()%100)==0) { 411 int nLook = model_>numberRows()+numberColumns; 412 number=0; 413 for (i=0;i<nLook;i++) { 414 if (infeas[i]) { 415 if (fabs(infeas[i])>COIN_INDEXED_TINY_ELEMENT) 416 index[number++]=i; 417 else 418 infeas[i]=0.0; 419 } 420 } 421 infeasible_>setNumElements(number); 422 } 423 if(model_>numberIterations()<model_>lastBadIteration()+200) { 424 // we can't really trust infeasibilities if there is dual error 425 double checkTolerance = 1.0e8; 426 if (!model_>factorization()>pivots()) 427 checkTolerance = 1.0e6; 428 if (model_>largestDualError()>checkTolerance) 429 tolerance *= model_>largestDualError()/checkTolerance; 430 // But cap 431 tolerance = min(1000.0,tolerance); 432 } 433 #ifdef CLP_DEBUG 434 if (model_>numberDualInfeasibilities()==1) 435 printf("** %g %g %g %x %x %d\n",tolerance,model_>dualTolerance(), 436 model_>largestDualError(),model_,model_>messageHandler(), 437 number); 438 #endif 439 // stop last one coming immediately 440 double saveOutInfeasibility=0.0; 441 if (sequenceOut>=0) { 442 saveOutInfeasibility = infeas[sequenceOut]; 443 infeas[sequenceOut]=0.0; 444 } 445 tolerance *= tolerance; // as we are using squares 446 447 int iPass; 448 // Setup two passes 449 int start[4]; 450 start[1]=number; 451 start[2]=0; 452 double dstart = ((double) number) * CoinDrand48(); 453 start[0]=(int) dstart; 454 start[3]=start[0]; 455 //double largestWeight=0.0; 456 //double smallestWeight=1.0e100; 457 for (iPass=0;iPass<2;iPass++) { 458 int end = start[2*iPass+1]; 459 if (switchType<5) { 460 for (i=start[2*iPass];i<end;i++) { 461 iSequence = index[i]; 462 double value = infeas[iSequence]; 463 if (value>tolerance) { 464 double weight = weights_[iSequence]; 465 //weight=1.0; 466 if (value>bestDj*weight) { 467 // check flagged variable and correct dj 468 if (!model_>flagged(iSequence)) { 469 bestDj=value/weight; 470 bestSequence = iSequence; 471 } else { 472 // just to make sure we don't exit before got something 473 numberWanted++; 474 } 475 } 476 } 477 numberWanted; 478 if (!numberWanted) 479 break; 480 } 481 } else { 482 // Dantzig 483 for (i=start[2*iPass];i<end;i++) { 484 iSequence = index[i]; 485 double value = infeas[iSequence]; 486 if (value>tolerance) { 487 if (value>bestDj) { 488 // check flagged variable and correct dj 489 if (!model_>flagged(iSequence)) { 490 bestDj=value; 491 bestSequence = iSequence; 492 } else { 493 // just to make sure we don't exit before got something 494 numberWanted++; 495 } 496 } 497 } 498 numberWanted; 499 if (!numberWanted) 500 break; 501 } 502 } 503 if (!numberWanted) 504 break; 505 } 506 if (sequenceOut>=0) { 507 infeas[sequenceOut]=saveOutInfeasibility; 508 } 509 /*if (model_>numberIterations()%100==0) 510 printf("%d best %g\n",bestSequence,bestDj);*/ 511 512 #ifdef CLP_DEBUG 513 if (bestSequence>=0) { 514 if (model_>getStatus(bestSequence)==ClpSimplex::atLowerBound) 515 assert(model_>reducedCost(bestSequence)<0.0); 516 if (model_>getStatus(bestSequence)==ClpSimplex::atUpperBound) 517 assert(model_>reducedCost(bestSequence)>0.0); 518 } 519 #endif 520 return bestSequence; 521 } 522 // Just update djs 523 void 524 ClpPrimalColumnSteepest::justDjs(CoinIndexedVector * updates, 525 CoinIndexedVector * spareRow1, 526 CoinIndexedVector * spareRow2, 527 CoinIndexedVector * spareColumn1, 528 CoinIndexedVector * spareColumn2) 529 { 162 530 int iSection,j; 163 int number; 531 int number=0; 532 int * index; 533 double * updateBy; 534 double * reducedCost; 535 double tolerance=model_>currentDualTolerance(); 536 // we can't really trust infeasibilities if there is dual error 537 // this coding has to mimic coding in checkDualSolution 538 double error = min(1.0e3,model_>largestDualError()); 539 // allow tolerance at least slightly bigger than standard 540 tolerance = tolerance + error; 541 int pivotRow = model_>pivotRow(); 542 double * infeas = infeasible_>denseVector(); 543 model_>factorization()>updateColumnTranspose(spareRow2,updates); 544 545 // put row of tableau in rowArray and columnArray 546 model_>clpMatrix()>transposeTimes(model_,1.0, 547 updates,spareColumn2,spareColumn1); 548 // normal 549 for (iSection=0;iSection<2;iSection++) { 550 551 reducedCost=model_>djRegion(iSection); 552 int addSequence; 553 554 if (!iSection) { 555 number = updates>getNumElements(); 556 index = updates>getIndices(); 557 updateBy = updates>denseVector(); 558 addSequence = model_>numberColumns(); 559 } else { 560 number = spareColumn1>getNumElements(); 561 index = spareColumn1>getIndices(); 562 updateBy = spareColumn1>denseVector(); 563 addSequence = 0; 564 } 565 566 for (j=0;j<number;j++) { 567 int iSequence = index[j]; 568 double value = reducedCost[iSequence]; 569 value = updateBy[iSequence]; 570 updateBy[iSequence]=0.0; 571 reducedCost[iSequence] = value; 572 ClpSimplex::Status status = model_>getStatus(iSequence+addSequence); 573 574 switch(status) { 575 576 case ClpSimplex::basic: 577 infeasible_>zero(iSequence+addSequence); 578 case ClpSimplex::isFixed: 579 break; 580 case ClpSimplex::isFree: 581 case ClpSimplex::superBasic: 582 if (fabs(value)>FREE_ACCEPT*tolerance) { 583 // we are going to bias towards free (but only if reasonable) 584 value *= FREE_BIAS; 585 // store square in list 586 if (infeas[iSequence+addSequence]) 587 infeas[iSequence+addSequence] = value*value; // already there 588 else 589 infeasible_>quickAdd(iSequence+addSequence,value*value); 590 } else { 591 infeasible_>zero(iSequence+addSequence); 592 } 593 break; 594 case ClpSimplex::atUpperBound: 595 if (value>tolerance) { 596 // store square in list 597 if (infeas[iSequence+addSequence]) 598 infeas[iSequence+addSequence] = value*value; // already there 599 else 600 infeasible_>quickAdd(iSequence+addSequence,value*value); 601 } else { 602 infeasible_>zero(iSequence+addSequence); 603 } 604 break; 605 case ClpSimplex::atLowerBound: 606 if (value<tolerance) { 607 // store square in list 608 if (infeas[iSequence+addSequence]) 609 infeas[iSequence+addSequence] = value*value; // already there 610 else 611 infeasible_>quickAdd(iSequence+addSequence,value*value); 612 } else { 613 infeasible_>zero(iSequence+addSequence); 614 } 615 } 616 } 617 } 618 updates>setNumElements(0); 619 spareColumn1>setNumElements(0); 620 if (pivotRow>=0) { 621 // make sure infeasibility on incoming is 0.0 622 int sequenceIn = model_>sequenceIn(); 623 infeasible_>zero(sequenceIn); 624 } 625 } 626 // Update djs, weights for Devex 627 void 628 ClpPrimalColumnSteepest::djsAndDevex(CoinIndexedVector * updates, 629 CoinIndexedVector * spareRow1, 630 CoinIndexedVector * spareRow2, 631 CoinIndexedVector * spareColumn1, 632 CoinIndexedVector * spareColumn2) 633 { 634 int j; 635 int number=0; 636 int * index; 637 double * updateBy; 638 double * reducedCost; 639 // dj could be very small (or even zero  take care) 640 double dj = model_>dualIn(); 641 double tolerance=model_>currentDualTolerance(); 642 // we can't really trust infeasibilities if there is dual error 643 // this coding has to mimic coding in checkDualSolution 644 double error = min(1.0e3,model_>largestDualError()); 645 // allow tolerance at least slightly bigger than standard 646 tolerance = tolerance + error; 647 // for weights update we use pivotSequence 648 // unset in case sub flip 649 assert (pivotSequence_>=0); 650 pivotSequence_=1; 651 double * infeas = infeasible_>denseVector(); 652 model_>factorization()>updateColumnTranspose(spareRow2,updates); 653 // and we can see if reference 654 double referenceIn=0.0; 655 int sequenceIn = model_>sequenceIn(); 656 if (mode_!=1&&reference(sequenceIn)) 657 referenceIn=1.0; 658 // save outgoing weight round update 659 double outgoingWeight=0.0; 660 int sequenceOut = model_>sequenceOut(); 661 if (sequenceOut>=0) 662 outgoingWeight=weights_[sequenceOut]; 663 664 // put row of tableau in rowArray and columnArray 665 model_>clpMatrix()>transposeTimes(model_,1.0, 666 updates,spareColumn2,spareColumn1); 667 // update weights 668 double * weight; 669 int numberColumns = model_>numberColumns(); 670 double scaleFactor = 1.0/dj; // as formula is with 1.0 671 // rows 672 reducedCost=model_>djRegion(0); 673 int addSequence=model_>numberColumns();; 674 675 number = updates>getNumElements(); 676 index = updates>getIndices(); 677 updateBy = updates>denseVector(); 678 weight = weights_+numberColumns; 679 // Devex 680 for (j=0;j<number;j++) { 681 double thisWeight; 682 double pivot; 683 double value3; 684 int iSequence = index[j]; 685 double value = reducedCost[iSequence]; 686 double value2 = updateBy[iSequence]; 687 updateBy[iSequence]=0.0; 688 value = value2; 689 reducedCost[iSequence] = value; 690 ClpSimplex::Status status = model_>getStatus(iSequence+addSequence); 691 692 switch(status) { 693 694 case ClpSimplex::basic: 695 infeasible_>zero(iSequence+addSequence); 696 case ClpSimplex::isFixed: 697 break; 698 case ClpSimplex::isFree: 699 case ClpSimplex::superBasic: 700 thisWeight = weight[iSequence]; 701 // row has 1 702 pivot = value2*scaleFactor; 703 value3 = pivot * pivot*devex_; 704 if (reference(iSequence+numberColumns)) 705 value3 += 1.0; 706 weight[iSequence] = max(0.99*thisWeight,value3); 707 if (fabs(value)>FREE_ACCEPT*tolerance) { 708 // we are going to bias towards free (but only if reasonable) 709 value *= FREE_BIAS; 710 // store square in list 711 if (infeas[iSequence+addSequence]) 712 infeas[iSequence+addSequence] = value*value; // already there 713 else 714 infeasible_>quickAdd(iSequence+addSequence,value*value); 715 } else { 716 infeasible_>zero(iSequence+addSequence); 717 } 718 break; 719 case ClpSimplex::atUpperBound: 720 thisWeight = weight[iSequence]; 721 // row has 1 722 pivot = value2*scaleFactor; 723 value3 = pivot * pivot*devex_; 724 if (reference(iSequence+numberColumns)) 725 value3 += 1.0; 726 weight[iSequence] = max(0.99*thisWeight,value3); 727 if (value>tolerance) { 728 // store square in list 729 if (infeas[iSequence+addSequence]) 730 infeas[iSequence+addSequence] = value*value; // already there 731 else 732 infeasible_>quickAdd(iSequence+addSequence,value*value); 733 } else { 734 infeasible_>zero(iSequence+addSequence); 735 } 736 break; 737 case ClpSimplex::atLowerBound: 738 thisWeight = weight[iSequence]; 739 // row has 1 740 pivot = value2*scaleFactor; 741 value3 = pivot * pivot*devex_; 742 if (reference(iSequence+numberColumns)) 743 value3 += 1.0; 744 weight[iSequence] = max(0.99*thisWeight,value3); 745 if (value<tolerance) { 746 // store square in list 747 if (infeas[iSequence+addSequence]) 748 infeas[iSequence+addSequence] = value*value; // already there 749 else 750 infeasible_>quickAdd(iSequence+addSequence,value*value); 751 } else { 752 infeasible_>zero(iSequence+addSequence); 753 } 754 } 755 } 756 757 // columns 758 weight = weights_; 759 760 scaleFactor = scaleFactor; 761 reducedCost=model_>djRegion(1); 762 number = spareColumn1>getNumElements(); 763 index = spareColumn1>getIndices(); 764 updateBy = spareColumn1>denseVector(); 765 766 // Devex 767 768 for (j=0;j<number;j++) { 769 double thisWeight; 770 double pivot; 771 double value3; 772 int iSequence = index[j]; 773 double value = reducedCost[iSequence]; 774 double value2 = updateBy[iSequence]; 775 value = value2; 776 updateBy[iSequence]=0.0; 777 reducedCost[iSequence] = value; 778 ClpSimplex::Status status = model_>getStatus(iSequence); 779 780 switch(status) { 781 782 case ClpSimplex::basic: 783 infeasible_>zero(iSequence); 784 case ClpSimplex::isFixed: 785 break; 786 case ClpSimplex::isFree: 787 case ClpSimplex::superBasic: 788 thisWeight = weight[iSequence]; 789 // row has 1 790 pivot = value2*scaleFactor; 791 value3 = pivot * pivot*devex_; 792 if (reference(iSequence)) 793 value3 += 1.0; 794 weight[iSequence] = max(0.99*thisWeight,value3); 795 if (fabs(value)>FREE_ACCEPT*tolerance) { 796 // we are going to bias towards free (but only if reasonable) 797 value *= FREE_BIAS; 798 // store square in list 799 if (infeas[iSequence]) 800 infeas[iSequence] = value*value; // already there 801 else 802 infeasible_>quickAdd(iSequence,value*value); 803 } else { 804 infeasible_>zero(iSequence); 805 } 806 break; 807 case ClpSimplex::atUpperBound: 808 thisWeight = weight[iSequence]; 809 // row has 1 810 pivot = value2*scaleFactor; 811 value3 = pivot * pivot*devex_; 812 if (reference(iSequence)) 813 value3 += 1.0; 814 weight[iSequence] = max(0.99*thisWeight,value3); 815 if (value>tolerance) { 816 // store square in list 817 if (infeas[iSequence]) 818 infeas[iSequence] = value*value; // already there 819 else 820 infeasible_>quickAdd(iSequence,value*value); 821 } else { 822 infeasible_>zero(iSequence); 823 } 824 break; 825 case ClpSimplex::atLowerBound: 826 thisWeight = weight[iSequence]; 827 // row has 1 828 pivot = value2*scaleFactor; 829 value3 = pivot * pivot*devex_; 830 if (reference(iSequence)) 831 value3 += 1.0; 832 weight[iSequence] = max(0.99*thisWeight,value3); 833 if (value<tolerance) { 834 // store square in list 835 if (infeas[iSequence]) 836 infeas[iSequence] = value*value; // already there 837 else 838 infeasible_>quickAdd(iSequence,value*value); 839 } else { 840 infeasible_>zero(iSequence); 841 } 842 } 843 } 844 // restore outgoing weight 845 if (sequenceOut>=0) 846 weights_[sequenceOut]=outgoingWeight; 847 // make sure infeasibility on incoming is 0.0 848 infeasible_>zero(sequenceIn); 849 spareRow2>setNumElements(0); 850 //#define SOME_DEBUG_1 851 #ifdef SOME_DEBUG_1 852 // check for accuracy 853 int iCheck=229; 854 //printf("weight for iCheck is %g\n",weights_[iCheck]); 855 int numberRows = model_>numberRows(); 856 //int numberColumns = model_>numberColumns(); 857 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 858 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 859 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 860 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 861 } 862 #endif 863 updates>setNumElements(0); 864 spareColumn1>setNumElements(0); 865 } 866 // Update djs, weights for Steepest 867 void 868 ClpPrimalColumnSteepest::djsAndSteepest(CoinIndexedVector * updates, 869 CoinIndexedVector * spareRow1, 870 CoinIndexedVector * spareRow2, 871 CoinIndexedVector * spareColumn1, 872 CoinIndexedVector * spareColumn2) 873 { 874 int j; 875 int number=0; 876 int * index; 877 double * updateBy; 878 double * reducedCost; 879 // dj could be very small (or even zero  take care) 880 double dj = model_>dualIn(); 881 double tolerance=model_>currentDualTolerance(); 882 // we can't really trust infeasibilities if there is dual error 883 // this coding has to mimic coding in checkDualSolution 884 double error = min(1.0e3,model_>largestDualError()); 885 // allow tolerance at least slightly bigger than standard 886 tolerance = tolerance + error; 887 // for weights update we use pivotSequence 888 // unset in case sub flip 889 assert (pivotSequence_>=0); 890 pivotSequence_=1; 891 double * infeas = infeasible_>denseVector(); 892 model_>factorization()>updateColumnTranspose(spareRow2,updates); 893 894 model_>factorization()>updateColumnTranspose(spareRow2, 895 alternateWeights_); 896 // and we can see if reference 897 double referenceIn=0.0; 898 int sequenceIn = model_>sequenceIn(); 899 if (mode_!=1&&reference(sequenceIn)) 900 referenceIn=1.0; 901 // save outgoing weight round update 902 double outgoingWeight=0.0; 903 int sequenceOut = model_>sequenceOut(); 904 if (sequenceOut>=0) 905 outgoingWeight=weights_[sequenceOut]; 906 907 // put row of tableau in rowArray and columnArray 908 model_>clpMatrix()>transposeTimes(model_,1.0, 909 updates,spareColumn2,spareColumn1); 910 // get subset which have nonzero tableau elements 911 // Luckily spareRow2 is long enough (rowArray_[3]) 912 model_>clpMatrix()>subsetTransposeTimes(model_,alternateWeights_, 913 spareColumn1, 914 spareRow2); 915 // update weights 916 double * weight; 917 double * other = alternateWeights_>denseVector(); 918 int numberColumns = model_>numberColumns(); 919 double scaleFactor = 1.0/dj; // as formula is with 1.0 920 // rows 921 reducedCost=model_>djRegion(0); 922 int addSequence=model_>numberColumns();; 923 924 number = updates>getNumElements(); 925 index = updates>getIndices(); 926 updateBy = updates>denseVector(); 927 weight = weights_+numberColumns; 928 929 for (j=0;j<number;j++) { 930 double thisWeight; 931 double pivot; 932 double modification; 933 double pivotSquared; 934 int iSequence = index[j]; 935 double value = reducedCost[iSequence]; 936 double value2 = updateBy[iSequence]; 937 value = value2; 938 reducedCost[iSequence] = value; 939 ClpSimplex::Status status = model_>getStatus(iSequence+addSequence); 940 updateBy[iSequence]=0.0; 941 942 switch(status) { 943 944 case ClpSimplex::basic: 945 infeasible_>zero(iSequence+addSequence); 946 case ClpSimplex::isFixed: 947 break; 948 case ClpSimplex::isFree: 949 case ClpSimplex::superBasic: 950 thisWeight = weight[iSequence]; 951 // row has 1 952 pivot = value2*scaleFactor; 953 modification = other[iSequence]; 954 pivotSquared = pivot * pivot; 955 956 thisWeight += pivotSquared * devex_ + pivot * modification; 957 if (thisWeight<TRY_NORM) { 958 if (mode_==1) { 959 // steepest 960 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 961 } else { 962 // exact 963 thisWeight = referenceIn*pivotSquared; 964 if (reference(iSequence+numberColumns)) 965 thisWeight += 1.0; 966 thisWeight = max(thisWeight,TRY_NORM); 967 } 968 } 969 weight[iSequence] = thisWeight; 970 if (fabs(value)>FREE_ACCEPT*tolerance) { 971 // we are going to bias towards free (but only if reasonable) 972 value *= FREE_BIAS; 973 // store square in list 974 if (infeas[iSequence+addSequence]) 975 infeas[iSequence+addSequence] = value*value; // already there 976 else 977 infeasible_>quickAdd(iSequence+addSequence,value*value); 978 } else { 979 infeasible_>zero(iSequence+addSequence); 980 } 981 break; 982 case ClpSimplex::atUpperBound: 983 thisWeight = weight[iSequence]; 984 // row has 1 985 pivot = value2*scaleFactor; 986 modification = other[iSequence]; 987 pivotSquared = pivot * pivot; 988 989 thisWeight += pivotSquared * devex_ + pivot * modification; 990 if (thisWeight<TRY_NORM) { 991 if (mode_==1) { 992 // steepest 993 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 994 } else { 995 // exact 996 thisWeight = referenceIn*pivotSquared; 997 if (reference(iSequence+numberColumns)) 998 thisWeight += 1.0; 999 thisWeight = max(thisWeight,TRY_NORM); 1000 } 1001 } 1002 weight[iSequence] = thisWeight; 1003 if (value>tolerance) { 1004 // store square in list 1005 if (infeas[iSequence+addSequence]) 1006 infeas[iSequence+addSequence] = value*value; // already there 1007 else 1008 infeasible_>quickAdd(iSequence+addSequence,value*value); 1009 } else { 1010 infeasible_>zero(iSequence+addSequence); 1011 } 1012 break; 1013 case ClpSimplex::atLowerBound: 1014 thisWeight = weight[iSequence]; 1015 // row has 1 1016 pivot = value2*scaleFactor; 1017 modification = other[iSequence]; 1018 pivotSquared = pivot * pivot; 1019 1020 thisWeight += pivotSquared * devex_ + pivot * modification; 1021 if (thisWeight<TRY_NORM) { 1022 if (mode_==1) { 1023 // steepest 1024 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1025 } else { 1026 // exact 1027 thisWeight = referenceIn*pivotSquared; 1028 if (reference(iSequence+numberColumns)) 1029 thisWeight += 1.0; 1030 thisWeight = max(thisWeight,TRY_NORM); 1031 } 1032 } 1033 weight[iSequence] = thisWeight; 1034 if (value<tolerance) { 1035 // store square in list 1036 if (infeas[iSequence+addSequence]) 1037 infeas[iSequence+addSequence] = value*value; // already there 1038 else 1039 infeasible_>quickAdd(iSequence+addSequence,value*value); 1040 } else { 1041 infeasible_>zero(iSequence+addSequence); 1042 } 1043 } 1044 } 1045 1046 // columns 1047 weight = weights_; 1048 1049 scaleFactor = scaleFactor; 1050 reducedCost=model_>djRegion(1); 1051 number = spareColumn1>getNumElements(); 1052 index = spareColumn1>getIndices(); 1053 updateBy = spareColumn1>denseVector(); 1054 double * updateBy2 = spareRow2>denseVector(); 1055 1056 for (j=0;j<number;j++) { 1057 double thisWeight; 1058 double pivot; 1059 double modification; 1060 double pivotSquared; 1061 int iSequence = index[j]; 1062 double value = reducedCost[iSequence]; 1063 double value2 = updateBy[iSequence]; 1064 updateBy[iSequence]=0.0; 1065 value = value2; 1066 reducedCost[iSequence] = value; 1067 ClpSimplex::Status status = model_>getStatus(iSequence); 1068 1069 switch(status) { 1070 1071 case ClpSimplex::basic: 1072 infeasible_>zero(iSequence); 1073 case ClpSimplex::isFixed: 1074 updateBy2[iSequence]=0.0; 1075 break; 1076 case ClpSimplex::isFree: 1077 case ClpSimplex::superBasic: 1078 thisWeight = weight[iSequence]; 1079 pivot = value2*scaleFactor; 1080 modification = updateBy2[iSequence]; 1081 updateBy2[iSequence]=0.0; 1082 pivotSquared = pivot * pivot; 1083 1084 thisWeight += pivotSquared * devex_ + pivot * modification; 1085 if (thisWeight<TRY_NORM) { 1086 if (mode_==1) { 1087 // steepest 1088 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1089 } else { 1090 // exact 1091 thisWeight = referenceIn*pivotSquared; 1092 if (reference(iSequence)) 1093 thisWeight += 1.0; 1094 thisWeight = max(thisWeight,TRY_NORM); 1095 } 1096 } 1097 weight[iSequence] = thisWeight; 1098 if (fabs(value)>FREE_ACCEPT*tolerance) { 1099 // we are going to bias towards free (but only if reasonable) 1100 value *= FREE_BIAS; 1101 // store square in list 1102 if (infeas[iSequence]) 1103 infeas[iSequence] = value*value; // already there 1104 else 1105 infeasible_>quickAdd(iSequence,value*value); 1106 } else { 1107 infeasible_>zero(iSequence); 1108 } 1109 break; 1110 case ClpSimplex::atUpperBound: 1111 thisWeight = weight[iSequence]; 1112 pivot = value2*scaleFactor; 1113 modification = updateBy2[iSequence]; 1114 updateBy2[iSequence]=0.0; 1115 pivotSquared = pivot * pivot; 1116 1117 thisWeight += pivotSquared * devex_ + pivot * modification; 1118 if (thisWeight<TRY_NORM) { 1119 if (mode_==1) { 1120 // steepest 1121 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1122 } else { 1123 // exact 1124 thisWeight = referenceIn*pivotSquared; 1125 if (reference(iSequence)) 1126 thisWeight += 1.0; 1127 thisWeight = max(thisWeight,TRY_NORM); 1128 } 1129 } 1130 weight[iSequence] = thisWeight; 1131 if (value>tolerance) { 1132 // store square in list 1133 if (infeas[iSequence]) 1134 infeas[iSequence] = value*value; // already there 1135 else 1136 infeasible_>quickAdd(iSequence,value*value); 1137 } else { 1138 infeasible_>zero(iSequence); 1139 } 1140 break; 1141 case ClpSimplex::atLowerBound: 1142 thisWeight = weight[iSequence]; 1143 pivot = value2*scaleFactor; 1144 modification = updateBy2[iSequence]; 1145 updateBy2[iSequence]=0.0; 1146 pivotSquared = pivot * pivot; 1147 1148 thisWeight += pivotSquared * devex_ + pivot * modification; 1149 if (thisWeight<TRY_NORM) { 1150 if (mode_==1) { 1151 // steepest 1152 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1153 } else { 1154 // exact 1155 thisWeight = referenceIn*pivotSquared; 1156 if (reference(iSequence)) 1157 thisWeight += 1.0; 1158 thisWeight = max(thisWeight,TRY_NORM); 1159 } 1160 } 1161 weight[iSequence] = thisWeight; 1162 if (value<tolerance) { 1163 // store square in list 1164 if (infeas[iSequence]) 1165 infeas[iSequence] = value*value; // already there 1166 else 1167 infeasible_>quickAdd(iSequence,value*value); 1168 } else { 1169 infeasible_>zero(iSequence); 1170 } 1171 } 1172 } 1173 // restore outgoing weight 1174 if (sequenceOut>=0) 1175 weights_[sequenceOut]=outgoingWeight; 1176 // make sure infeasibility on incoming is 0.0 1177 infeasible_>zero(sequenceIn); 1178 alternateWeights_>clear(); 1179 spareRow2>setNumElements(0); 1180 //#define SOME_DEBUG_1 1181 #ifdef SOME_DEBUG_1 1182 // check for accuracy 1183 int iCheck=229; 1184 //printf("weight for iCheck is %g\n",weights_[iCheck]); 1185 int numberRows = model_>numberRows(); 1186 //int numberColumns = model_>numberColumns(); 1187 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 1188 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 1189 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 1190 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 1191 } 1192 #endif 1193 updates>setNumElements(0); 1194 spareColumn1>setNumElements(0); 1195 } 1196 // Update djs, weights for Devex 1197 void 1198 ClpPrimalColumnSteepest::djsAndDevex2(CoinIndexedVector * updates, 1199 CoinIndexedVector * spareRow1, 1200 CoinIndexedVector * spareRow2, 1201 CoinIndexedVector * spareColumn1, 1202 CoinIndexedVector * spareColumn2) 1203 { 1204 int iSection,j; 1205 int number=0; 164 1206 int * index; 165 1207 double * updateBy; … … 174 1216 tolerance = tolerance + error; 175 1217 int pivotRow = model_>pivotRow(); 176 // Local copy of mode so can decide what to do 177 int switchType = mode_; 178 if (mode_==4) { 179 if (numberSwitched_) { 180 switchType=3; 1218 double * infeas = infeasible_>denseVector(); 1219 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1220 1221 // put row of tableau in rowArray and columnArray 1222 model_>clpMatrix()>transposeTimes(model_,1.0, 1223 updates,spareColumn2,spareColumn1); 1224 // normal 1225 for (iSection=0;iSection<2;iSection++) { 1226 1227 reducedCost=model_>djRegion(iSection); 1228 int addSequence; 1229 1230 if (!iSection) { 1231 number = updates>getNumElements(); 1232 index = updates>getIndices(); 1233 updateBy = updates>denseVector(); 1234 addSequence = model_>numberColumns(); 181 1235 } else { 182 // See if to switch 183 int numberElements = model_>factorization()>numberElements(); 184 int numberRows = model_>numberRows(); 185 // ratio is done on number of columns here 186 int numberColumns = model_>numberColumns(); 187 double ratio = (double) numberElements/(double) numberColumns; 188 //double ratio = (double) numberElements/(double) model_>clpMatrix()>getNumElements(); 189 int numberWanted; 190 bool doPartial=true; 191 int numberTotal = numberColumns+model_>numberRows(); 192 if (ratio<0.1) { 193 numberWanted = max(100,numberTotal/200); 194 } else if (ratio<1.0) { 195 numberWanted = max(100,numberTotal/20); 196 } else if (ratio<2.0) { 197 numberWanted = max(200,numberTotal/10); 1236 number = spareColumn1>getNumElements(); 1237 index = spareColumn1>getIndices(); 1238 updateBy = spareColumn1>denseVector(); 1239 addSequence = 0; 1240 } 1241 1242 for (j=0;j<number;j++) { 1243 int iSequence = index[j]; 1244 double value = reducedCost[iSequence]; 1245 value = updateBy[iSequence]; 1246 reducedCost[iSequence] = value; 1247 ClpSimplex::Status status = model_>getStatus(iSequence+addSequence); 1248 1249 switch(status) { 1250 1251 case ClpSimplex::basic: 1252 infeasible_>zero(iSequence+addSequence); 1253 case ClpSimplex::isFixed: 1254 break; 1255 case ClpSimplex::isFree: 1256 case ClpSimplex::superBasic: 1257 if (fabs(value)>FREE_ACCEPT*tolerance) { 1258 // we are going to bias towards free (but only if reasonable) 1259 value *= FREE_BIAS; 1260 // store square in list 1261 if (infeas[iSequence+addSequence]) 1262 infeas[iSequence+addSequence] = value*value; // already there 1263 else 1264 infeasible_>quickAdd(iSequence+addSequence,value*value); 1265 } else { 1266 infeasible_>zero(iSequence+addSequence); 1267 } 1268 break; 1269 case ClpSimplex::atUpperBound: 1270 if (value>tolerance) { 1271 // store square in list 1272 if (infeas[iSequence+addSequence]) 1273 infeas[iSequence+addSequence] = value*value; // already there 1274 else 1275 infeasible_>quickAdd(iSequence+addSequence,value*value); 1276 } else { 1277 infeasible_>zero(iSequence+addSequence); 1278 } 1279 break; 1280 case ClpSimplex::atLowerBound: 1281 if (value<tolerance) { 1282 // store square in list 1283 if (infeas[iSequence+addSequence]) 1284 infeas[iSequence+addSequence] = value*value; // already there 1285 else 1286 infeasible_>quickAdd(iSequence+addSequence,value*value); 1287 } else { 1288 infeasible_>zero(iSequence+addSequence); 1289 } 1290 } 1291 } 1292 } 1293 // we can zero out as will have to get pivot row 1294 updates>clear(); 1295 spareColumn1>clear(); 1296 // make sure infeasibility on incoming is 0.0 1297 int sequenceIn = model_>sequenceIn(); 1298 infeasible_>zero(sequenceIn); 1299 // for weights update we use pivotSequence 1300 assert (pivotSequence_>=0); 1301 pivotRow = pivotSequence_; 1302 // unset in case sub flip 1303 pivotSequence_=1; 1304 // make sure infeasibility on incoming is 0.0 1305 const int * pivotVariable = model_>pivotVariable(); 1306 sequenceIn = pivotVariable[pivotRow]; 1307 infeasible_>zero(sequenceIn); 1308 // and we can see if reference 1309 double referenceIn=0.0; 1310 if (mode_!=1&&reference(sequenceIn)) 1311 referenceIn=1.0; 1312 // save outgoing weight round update 1313 double outgoingWeight=0.0; 1314 int sequenceOut = model_>sequenceOut(); 1315 if (sequenceOut>=0) 1316 outgoingWeight=weights_[sequenceOut]; 1317 // update weights 1318 updates>setNumElements(0); 1319 spareColumn1>setNumElements(0); 1320 // might as well set dj to 1 1321 dj = 1.0; 1322 updates>insert(pivotRow,dj); 1323 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1324 // put row of tableau in rowArray and columnArray 1325 model_>clpMatrix()>transposeTimes(model_,1.0, 1326 updates,spareColumn2,spareColumn1); 1327 double * weight; 1328 int numberColumns = model_>numberColumns(); 1329 double scaleFactor = 1.0/dj; // as formula is with 1.0 1330 // rows 1331 number = updates>getNumElements(); 1332 index = updates>getIndices(); 1333 updateBy = updates>denseVector(); 1334 weight = weights_+numberColumns; 1335 1336 assert (devex_>0.0); 1337 for (j=0;j<number;j++) { 1338 int iSequence = index[j]; 1339 double thisWeight = weight[iSequence]; 1340 // row has 1 1341 double pivot = updateBy[iSequence]*scaleFactor; 1342 updateBy[iSequence]=0.0; 1343 double value = pivot * pivot*devex_; 1344 if (reference(iSequence+numberColumns)) 1345 value += 1.0; 1346 weight[iSequence] = max(0.99*thisWeight,value); 1347 } 1348 1349 // columns 1350 weight = weights_; 1351 1352 scaleFactor = scaleFactor; 1353 1354 number = spareColumn1>getNumElements(); 1355 index = spareColumn1>getIndices(); 1356 updateBy = spareColumn1>denseVector(); 1357 for (j=0;j<number;j++) { 1358 int iSequence = index[j]; 1359 double thisWeight = weight[iSequence]; 1360 // row has 1 1361 double pivot = updateBy[iSequence]*scaleFactor; 1362 updateBy[iSequence]=0.0; 1363 double value = pivot * pivot*devex_; 1364 if (reference(iSequence)) 1365 value += 1.0; 1366 weight[iSequence] = max(0.99*thisWeight,value); 1367 } 1368 // restore outgoing weight 1369 if (sequenceOut>=0) 1370 weights_[sequenceOut]=outgoingWeight; 1371 spareColumn2>setNumElements(0); 1372 //#define SOME_DEBUG_1 1373 #ifdef SOME_DEBUG_1 1374 // check for accuracy 1375 int iCheck=229; 1376 //printf("weight for iCheck is %g\n",weights_[iCheck]); 1377 int numberRows = model_>numberRows(); 1378 //int numberColumns = model_>numberColumns(); 1379 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 1380 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 1381 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 1382 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 1383 } 1384 #endif 1385 updates>setNumElements(0); 1386 spareColumn1>setNumElements(0); 1387 } 1388 // Update djs, weights for Steepest 1389 void 1390 ClpPrimalColumnSteepest::djsAndSteepest2(CoinIndexedVector * updates, 1391 CoinIndexedVector * spareRow1, 1392 CoinIndexedVector * spareRow2, 1393 CoinIndexedVector * spareColumn1, 1394 CoinIndexedVector * spareColumn2) 1395 { 1396 int iSection,j; 1397 int number=0; 1398 int * index; 1399 double * updateBy; 1400 double * reducedCost; 1401 // dj could be very small (or even zero  take care) 1402 double dj = model_>dualIn(); 1403 double tolerance=model_>currentDualTolerance(); 1404 // we can't really trust infeasibilities if there is dual error 1405 // this coding has to mimic coding in checkDualSolution 1406 double error = min(1.0e3,model_>largestDualError()); 1407 // allow tolerance at least slightly bigger than standard 1408 tolerance = tolerance + error; 1409 int pivotRow = model_>pivotRow(); 1410 double * infeas = infeasible_>denseVector(); 1411 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1412 1413 // put row of tableau in rowArray and columnArray 1414 model_>clpMatrix()>transposeTimes(model_,1.0, 1415 updates,spareColumn2,spareColumn1); 1416 // normal 1417 for (iSection=0;iSection<2;iSection++) { 1418 1419 reducedCost=model_>djRegion(iSection); 1420 int addSequence; 1421 1422 if (!iSection) { 1423 number = updates>getNumElements(); 1424 index = updates>getIndices(); 1425 updateBy = updates>denseVector(); 1426 addSequence = model_>numberColumns(); 1427 } else { 1428 number = spareColumn1>getNumElements(); 1429 index = spareColumn1>getIndices(); 1430 updateBy = spareColumn1>denseVector(); 1431 addSequence = 0; 1432 } 1433 1434 for (j=0;j<number;j++) { 1435 int iSequence = index[j]; 1436 double value = reducedCost[iSequence]; 1437 value = updateBy[iSequence]; 1438 reducedCost[iSequence] = value; 1439 ClpSimplex::Status status = model_>getStatus(iSequence+addSequence); 1440 1441 switch(status) { 1442 1443 case ClpSimplex::basic: 1444 infeasible_>zero(iSequence+addSequence); 1445 case ClpSimplex::isFixed: 1446 break; 1447 case ClpSimplex::isFree: 1448 case ClpSimplex::superBasic: 1449 if (fabs(value)>FREE_ACCEPT*tolerance) { 1450 // we are going to bias towards free (but only if reasonable) 1451 value *= FREE_BIAS; 1452 // store square in list 1453 if (infeas[iSequence+addSequence]) 1454 infeas[iSequence+addSequence] = value*value; // already there 1455 else 1456 infeasible_>quickAdd(iSequence+addSequence,value*value); 1457 } else { 1458 infeasible_>zero(iSequence+addSequence); 1459 } 1460 break; 1461 case ClpSimplex::atUpperBound: 1462 if (value>tolerance) { 1463 // store square in list 1464 if (infeas[iSequence+addSequence]) 1465 infeas[iSequence+addSequence] = value*value; // already there 1466 else 1467 infeasible_>quickAdd(iSequence+addSequence,value*value); 1468 } else { 1469 infeasible_>zero(iSequence+addSequence); 1470 } 1471 break; 1472 case ClpSimplex::atLowerBound: 1473 if (value<tolerance) { 1474 // store square in list 1475 if (infeas[iSequence+addSequence]) 1476 infeas[iSequence+addSequence] = value*value; // already there 1477 else 1478 infeasible_>quickAdd(iSequence+addSequence,value*value); 1479 } else { 1480 infeasible_>zero(iSequence+addSequence); 1481 } 1482 } 1483 } 1484 } 1485 // we can zero out as will have to get pivot row 1486 // ***** move 1487 updates>clear(); 1488 spareColumn1>clear(); 1489 if (pivotRow>=0) { 1490 // make sure infeasibility on incoming is 0.0 1491 int sequenceIn = model_>sequenceIn(); 1492 infeasible_>zero(sequenceIn); 1493 } 1494 // for weights update we use pivotSequence 1495 pivotRow = pivotSequence_; 1496 // unset in case sub flip 1497 pivotSequence_=1; 1498 if (pivotRow>=0) { 1499 // make sure infeasibility on incoming is 0.0 1500 const int * pivotVariable = model_>pivotVariable(); 1501 int sequenceIn = pivotVariable[pivotRow]; 1502 infeasible_>zero(sequenceIn); 1503 // and we can see if reference 1504 double referenceIn=0.0; 1505 if (mode_!=1&&reference(sequenceIn)) 1506 referenceIn=1.0; 1507 // save outgoing weight round update 1508 double outgoingWeight=0.0; 1509 int sequenceOut = model_>sequenceOut(); 1510 if (sequenceOut>=0) 1511 outgoingWeight=weights_[sequenceOut]; 1512 // update weights 1513 updates>setNumElements(0); 1514 spareColumn1>setNumElements(0); 1515 // might as well set dj to 1 1516 dj = 1.0; 1517 updates>insert(pivotRow,dj); 1518 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1519 // put row of tableau in rowArray and columnArray 1520 model_>clpMatrix()>transposeTimes(model_,1.0, 1521 updates,spareColumn2,spareColumn1); 1522 double * weight; 1523 double * other = alternateWeights_>denseVector(); 1524 int numberColumns = model_>numberColumns(); 1525 double scaleFactor = 1.0/dj; // as formula is with 1.0 1526 // rows 1527 number = updates>getNumElements(); 1528 index = updates>getIndices(); 1529 updateBy = updates>denseVector(); 1530 weight = weights_+numberColumns; 1531 1532 if (mode_<4numberSwitched_>1) { 1533 // Exact 1534 // now update weight update array 1535 model_>factorization()>updateColumnTranspose(spareRow2, 1536 alternateWeights_); 1537 for (j=0;j<number;j++) { 1538 int iSequence = index[j]; 1539 double thisWeight = weight[iSequence]; 1540 // row has 1 1541 double pivot = updateBy[iSequence]*scaleFactor; 1542 updateBy[iSequence]=0.0; 1543 double modification = other[iSequence]; 1544 double pivotSquared = pivot * pivot; 1545 1546 thisWeight += pivotSquared * devex_ + pivot * modification; 1547 if (thisWeight<TRY_NORM) { 1548 if (mode_==1) { 1549 // steepest 1550 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1551 } else { 1552 // exact 1553 thisWeight = referenceIn*pivotSquared; 1554 if (reference(iSequence+numberColumns)) 1555 thisWeight += 1.0; 1556 thisWeight = max(thisWeight,TRY_NORM); 1557 } 1558 } 1559 weight[iSequence] = thisWeight; 1560 } 1561 } else if (mode_==4) { 1562 // Devex 1563 assert (devex_>0.0); 1564 for (j=0;j<number;j++) { 1565 int iSequence = index[j]; 1566 double thisWeight = weight[iSequence]; 1567 // row has 1 1568 double pivot = updateBy[iSequence]*scaleFactor; 1569 updateBy[iSequence]=0.0; 1570 double value = pivot * pivot*devex_; 1571 if (reference(iSequence+numberColumns)) 1572 value += 1.0; 1573 weight[iSequence] = max(0.99*thisWeight,value); 1574 } 1575 } 1576 1577 // columns 1578 weight = weights_; 1579 1580 scaleFactor = scaleFactor; 1581 1582 number = spareColumn1>getNumElements(); 1583 index = spareColumn1>getIndices(); 1584 updateBy = spareColumn1>denseVector(); 1585 if (mode_<4numberSwitched_>1) { 1586 // Exact 1587 // get subset which have nonzero tableau elements 1588 model_>clpMatrix()>subsetTransposeTimes(model_,alternateWeights_, 1589 spareColumn1, 1590 spareColumn2); 1591 double * updateBy2 = spareColumn2>denseVector(); 1592 for (j=0;j<number;j++) { 1593 int iSequence = index[j]; 1594 double thisWeight = weight[iSequence]; 1595 double pivot = updateBy[iSequence]*scaleFactor; 1596 updateBy[iSequence]=0.0; 1597 double modification = updateBy2[iSequence]; 1598 updateBy2[iSequence]=0.0; 1599 double pivotSquared = pivot * pivot; 1600 1601 thisWeight += pivotSquared * devex_ + pivot * modification; 1602 if (thisWeight<TRY_NORM) { 1603 if (mode_==1) { 1604 // steepest 1605 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1606 } else { 1607 // exact 1608 thisWeight = referenceIn*pivotSquared; 1609 if (reference(iSequence)) 1610 thisWeight += 1.0; 1611 thisWeight = max(thisWeight,TRY_NORM); 1612 } 1613 } 1614 weight[iSequence] = thisWeight; 1615 } 1616 } else if (mode_==4) { 1617 // Devex 1618 for (j=0;j<number;j++) { 1619 int iSequence = index[j]; 1620 double thisWeight = weight[iSequence]; 1621 // row has 1 1622 double pivot = updateBy[iSequence]*scaleFactor; 1623 updateBy[iSequence]=0.0; 1624 double value = pivot * pivot*devex_; 1625 if (reference(iSequence)) 1626 value += 1.0; 1627 weight[iSequence] = max(0.99*thisWeight,value); 1628 } 1629 } 1630 // restore outgoing weight 1631 if (sequenceOut>=0) 1632 weights_[sequenceOut]=outgoingWeight; 1633 alternateWeights_>clear(); 1634 spareColumn2>setNumElements(0); 1635 //#define SOME_DEBUG_1 1636 #ifdef SOME_DEBUG_1 1637 // check for accuracy 1638 int iCheck=229; 1639 //printf("weight for iCheck is %g\n",weights_[iCheck]); 1640 int numberRows = model_>numberRows(); 1641 //int numberColumns = model_>numberColumns(); 1642 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 1643 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 1644 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 1645 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 1646 } 1647 #endif 1648 } 1649 updates>setNumElements(0); 1650 spareColumn1>setNumElements(0); 1651 } 1652 // Update weights for Devex 1653 void 1654 ClpPrimalColumnSteepest::justDevex(CoinIndexedVector * updates, 1655 CoinIndexedVector * spareRow1, 1656 CoinIndexedVector * spareRow2, 1657 CoinIndexedVector * spareColumn1, 1658 CoinIndexedVector * spareColumn2) 1659 { 1660 int j; 1661 int number=0; 1662 int * index; 1663 double * updateBy; 1664 // dj could be very small (or even zero  take care) 1665 double dj = model_>dualIn(); 1666 double tolerance=model_>currentDualTolerance(); 1667 // we can't really trust infeasibilities if there is dual error 1668 // this coding has to mimic coding in checkDualSolution 1669 double error = min(1.0e3,model_>largestDualError()); 1670 // allow tolerance at least slightly bigger than standard 1671 tolerance = tolerance + error; 1672 int pivotRow = model_>pivotRow(); 1673 // for weights update we use pivotSequence 1674 pivotRow = pivotSequence_; 1675 assert (pivotRow>=0); 1676 // make sure infeasibility on incoming is 0.0 1677 const int * pivotVariable = model_>pivotVariable(); 1678 int sequenceIn = pivotVariable[pivotRow]; 1679 infeasible_>zero(sequenceIn); 1680 // and we can see if reference 1681 double referenceIn=0.0; 1682 if (mode_!=1&&reference(sequenceIn)) 1683 referenceIn=1.0; 1684 // save outgoing weight round update 1685 double outgoingWeight=0.0; 1686 int sequenceOut = model_>sequenceOut(); 1687 if (sequenceOut>=0) 1688 outgoingWeight=weights_[sequenceOut]; 1689 assert (!updates>getNumElements()); 1690 assert (!spareColumn1>getNumElements()); 1691 // unset in case sub flip 1692 pivotSequence_=1; 1693 // might as well set dj to 1 1694 dj = 1.0; 1695 updates>insert(pivotRow,dj); 1696 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1697 // put row of tableau in rowArray and columnArray 1698 model_>clpMatrix()>transposeTimes(model_,1.0, 1699 updates,spareColumn2,spareColumn1); 1700 double * weight; 1701 int numberColumns = model_>numberColumns(); 1702 double scaleFactor = 1.0/dj; // as formula is with 1.0 1703 // rows 1704 number = updates>getNumElements(); 1705 index = updates>getIndices(); 1706 updateBy = updates>denseVector(); 1707 weight = weights_+numberColumns; 1708 1709 // Devex 1710 assert (devex_>0.0); 1711 for (j=0;j<number;j++) { 1712 int iSequence = index[j]; 1713 double thisWeight = weight[iSequence]; 1714 // row has 1 1715 double pivot = updateBy[iSequence]*scaleFactor; 1716 updateBy[iSequence]=0.0; 1717 double value = pivot * pivot*devex_; 1718 if (reference(iSequence+numberColumns)) 1719 value += 1.0; 1720 weight[iSequence] = max(0.99*thisWeight,value); 1721 } 1722 1723 // columns 1724 weight = weights_; 1725 1726 scaleFactor = scaleFactor; 1727 1728 number = spareColumn1>getNumElements(); 1729 index = spareColumn1>getIndices(); 1730 updateBy = spareColumn1>denseVector(); 1731 // Devex 1732 for (j=0;j<number;j++) { 1733 int iSequence = index[j]; 1734 double thisWeight = weight[iSequence]; 1735 // row has 1 1736 double pivot = updateBy[iSequence]*scaleFactor; 1737 updateBy[iSequence]=0.0; 1738 double value = pivot * pivot*devex_; 1739 if (reference(iSequence)) 1740 value += 1.0; 1741 weight[iSequence] = max(0.99*thisWeight,value); 1742 } 1743 // restore outgoing weight 1744 if (sequenceOut>=0) 1745 weights_[sequenceOut]=outgoingWeight; 1746 spareColumn2>setNumElements(0); 1747 //#define SOME_DEBUG_1 1748 #ifdef SOME_DEBUG_1 1749 // check for accuracy 1750 int iCheck=229; 1751 //printf("weight for iCheck is %g\n",weights_[iCheck]); 1752 int numberRows = model_>numberRows(); 1753 //int numberColumns = model_>numberColumns(); 1754 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 1755 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 1756 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 1757 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 1758 } 1759 #endif 1760 updates>setNumElements(0); 1761 spareColumn1>setNumElements(0); 1762 } 1763 // Update weights for Steepest 1764 void 1765 ClpPrimalColumnSteepest::justSteepest(CoinIndexedVector * updates, 1766 CoinIndexedVector * spareRow1, 1767 CoinIndexedVector * spareRow2, 1768 CoinIndexedVector * spareColumn1, 1769 CoinIndexedVector * spareColumn2) 1770 { 1771 int j; 1772 int number=0; 1773 int * index; 1774 double * updateBy; 1775 // dj could be very small (or even zero  take care) 1776 double dj = model_>dualIn(); 1777 double tolerance=model_>currentDualTolerance(); 1778 // we can't really trust infeasibilities if there is dual error 1779 // this coding has to mimic coding in checkDualSolution 1780 double error = min(1.0e3,model_>largestDualError()); 1781 // allow tolerance at least slightly bigger than standard 1782 tolerance = tolerance + error; 1783 int pivotRow = model_>pivotRow(); 1784 // for weights update we use pivotSequence 1785 pivotRow = pivotSequence_; 1786 // unset in case sub flip 1787 pivotSequence_=1; 1788 assert (pivotRow>=0); 1789 // make sure infeasibility on incoming is 0.0 1790 const int * pivotVariable = model_>pivotVariable(); 1791 int sequenceIn = pivotVariable[pivotRow]; 1792 infeasible_>zero(sequenceIn); 1793 // and we can see if reference 1794 double referenceIn=0.0; 1795 if (mode_!=1&&reference(sequenceIn)) 1796 referenceIn=1.0; 1797 // save outgoing weight round update 1798 double outgoingWeight=0.0; 1799 int sequenceOut = model_>sequenceOut(); 1800 if (sequenceOut>=0) 1801 outgoingWeight=weights_[sequenceOut]; 1802 assert (!updates>getNumElements()); 1803 assert (!spareColumn1>getNumElements()); 1804 // update weights 1805 //updates>setNumElements(0); 1806 //spareColumn1>setNumElements(0); 1807 // might as well set dj to 1 1808 dj = 1.0; 1809 updates>insert(pivotRow,dj); 1810 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1811 // put row of tableau in rowArray and columnArray 1812 model_>clpMatrix()>transposeTimes(model_,1.0, 1813 updates,spareColumn2,spareColumn1); 1814 double * weight; 1815 double * other = alternateWeights_>denseVector(); 1816 int numberColumns = model_>numberColumns(); 1817 double scaleFactor = 1.0/dj; // as formula is with 1.0 1818 // rows 1819 number = updates>getNumElements(); 1820 index = updates>getIndices(); 1821 updateBy = updates>denseVector(); 1822 weight = weights_+numberColumns; 1823 1824 // Exact 1825 // now update weight update array 1826 model_>factorization()>updateColumnTranspose(spareRow2, 1827 alternateWeights_); 1828 for (j=0;j<number;j++) { 1829 int iSequence = index[j]; 1830 double thisWeight = weight[iSequence]; 1831 // row has 1 1832 double pivot = updateBy[iSequence]*scaleFactor; 1833 updateBy[iSequence]=0.0; 1834 double modification = other[iSequence]; 1835 double pivotSquared = pivot * pivot; 1836 1837 thisWeight += pivotSquared * devex_ + pivot * modification; 1838 if (thisWeight<TRY_NORM) { 1839 if (mode_==1) { 1840 // steepest 1841 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 198 1842 } else { 199 doPartial=false; // switch off 200 numberWanted=0; // just for compliler message 201 } 202 if (doPartial) { 203 if (model_>numberIterations()%100==0) 204 printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted); 205 // No do partial 206 // Always do slacks 207 int anyUpdates; 208 double * dual = model_>dualRowSolution(); 209 if (updates>getNumElements()) { 210 // would have to have two goes for devex, three for steepest 211 anyUpdates=2; 212 // add in pivot contribution 213 if (pivotRow>=0) 214 updates>add(pivotRow,dj); 215 } else if (pivotRow>=0) { 216 if (fabs(dj)>1.0e15) { 217 // some dj 218 updates>insert(pivotRow,dj); 219 if (fabs(dj)>1.0e6) { 220 // reasonable size 221 anyUpdates=1; 222 } else { 223 // too small 224 anyUpdates=2; 225 } 226 } else { 227 // zero dj 228 anyUpdates=1; 229 } 230 } else if (pivotSequence_>=0){ 231 // just after refactorization 232 anyUpdates=1; 233 } else { 234 // sub flip  nothing to do 235 anyUpdates=0; 236 } 237 238 // For now (as we can always switch off) 239 assert (!model_>nonLinearCost()>lookBothWays()); 240 double * slackDj =model_>djRegion(0); 241 double * dj =model_>djRegion(1); 242 if (anyUpdates>0) { 243 model_>factorization()>updateColumnTranspose(spareRow2,updates); 1843 // exact 1844 thisWeight = referenceIn*pivotSquared; 1845 if (reference(iSequence+numberColumns)) 1846 thisWeight += 1.0; 1847 thisWeight = max(thisWeight,TRY_NORM); 1848 } 1849 } 1850 weight[iSequence] = thisWeight; 1851 } 1852 1853 // columns 1854 weight = weights_; 1855 1856 scaleFactor = scaleFactor; 1857 1858 number = spareColumn1>getNumElements(); 1859 index = spareColumn1>getIndices(); 1860 updateBy = spareColumn1>denseVector(); 1861 // Exact 1862 // get subset which have nonzero tableau elements 1863 model_>clpMatrix()>subsetTransposeTimes(model_,alternateWeights_, 1864 spareColumn1, 1865 spareColumn2); 1866 double * updateBy2 = spareColumn2>denseVector(); 1867 for (j=0;j<number;j++) { 1868 int iSequence = index[j]; 1869 double thisWeight = weight[iSequence]; 1870 double pivot = updateBy[iSequence]*scaleFactor; 1871 updateBy[iSequence]=0.0; 1872 double modification = updateBy2[iSequence]; 1873 updateBy2[iSequence]=0.0; 1874 double pivotSquared = pivot * pivot; 244 1875 245 number = updates>getNumElements(); 246 index = updates>getIndices(); 247 updateBy = updates>denseVector(); 248 249 250 for (j=0;j<number;j++) { 251 int iSequence = index[j]; 252 double value = updateBy[iSequence]; 253 slackDj[iSequence] = value; 254 dual[iSequence] = value; 255 updateBy[iSequence]=0.0; 256 } 257 updates>setNumElements(0); 258 } 259 260 int iPass; 261 // Setup two passes AND treat slacks differently 262 int start[2][4]; 263 // slacks 264 start[0][1]=numberRows; 265 start[0][2]=0; 266 double dstart = ((double) numberRows) * CoinDrand48(); 267 start[0][0]=(int) dstart; 268 start[0][3]=start[0][0]; 269 for (iPass=0;iPass<4;iPass++) 270 start[0][iPass] += numberColumns; 271 start[1][1]=numberColumns; 272 start[1][2]=0; 273 dstart = ((double) numberColumns) * CoinDrand48(); 274 start[1][0]=(int) dstart; 275 start[1][3]=start[1][0]; 276 int bestSequence=1; 277 double bestDj=tolerance; 278 // We are going to fake it so pi looks sparse 279 double * saveDense = updates>denseVector(); 280 updates>setDenseVector(dual); 281 int kChunk = max(2000,numberWanted); 282 for (iPass=0;iPass<2;iPass++) { 283 // Slacks 284 int end = start[0][2*iPass+1]; 285 for (int iSequence=start[0][2*iPass];iSequence<end;iSequence++) { 286 double value; 287 ClpSimplex::Status status = model_>getStatus(iSequence); 288 289 switch(status) { 290 291 case ClpSimplex::basic: 292 case ClpSimplex::isFixed: 293 break; 294 case ClpSimplex::isFree: 295 case ClpSimplex::superBasic: 296 value = dj[iSequence]; 297 if (fabs(value)>FREE_ACCEPT*tolerance) { 298 numberWanted; 299 // we are going to bias towards free (but only if reasonable) 300 value = fabs(value*FREE_BIAS); 301 if (value>bestDj) { 302 // check flagged variable and correct dj 303 if (!model_>flagged(iSequence)) { 304 bestDj=value; 305 bestSequence = iSequence; 306 } else { 307 // just to make sure we don't exit before got something 308 numberWanted++; 309 } 310 } 311 } 312 break; 313 case ClpSimplex::atUpperBound: 314 value = dj[iSequence]; 315 if (value>tolerance) { 316 numberWanted; 317 if (value>bestDj) { 318 // check flagged variable and correct dj 319 if (!model_>flagged(iSequence)) { 320 bestDj=value; 321 bestSequence = iSequence; 322 } else { 323 // just to make sure we don't exit before got something 324 numberWanted++; 325 } 326 } 327 } 328 break; 329 case ClpSimplex::atLowerBound: 330 value = dj[iSequence]; 331 if (value<tolerance) { 332 numberWanted; 333 if (value>bestDj) { 334 // check flagged variable and correct dj 335 if (!model_>flagged(iSequence)) { 336 bestDj=value; 337 bestSequence = iSequence; 338 } else { 339 // just to make sure we don't exit before got something 340 numberWanted++; 341 } 342 } 343 } 344 } 345 if (!numberWanted) 346 break; 347 } 348 if (!numberWanted) 349 break; 350 // For Columns break up into manageable chunks 351 end = start[1][2*iPass+1]; 352 index = spareColumn1>getIndices(); 353 updateBy = spareColumn1>denseVector(); 354 for (int jSequence=start[1][2*iPass];jSequence<end;jSequence+=kChunk) { 355 int endThis = min(end,jSequence+kChunk); 356 number=0; 357 for (int iSequence=jSequence;iSequence<endThis;iSequence++) { 358 ClpSimplex::Status status = model_>getStatus(iSequence); 359 360 if(status!=ClpSimplex::basic&&status!=ClpSimplex::isFixed) { 361 index[number++]=iSequence; 362 } 363 } 364 spareColumn1>setNumElements(number); 365 // get subset which have nonzero djs 366 // updates is actually slack djs 367 model_>clpMatrix()>subsetTransposeTimes(model_,updates, 368 spareColumn1, 369 spareColumn2); 370 double * updateBy2 = spareColumn2>denseVector(); 371 const double * cost = model_>costRegion(); 372 for (int j=0;j<number;j++) { 373 int iSequence = index[j]; 374 double value = cost[iSequence]updateBy2[iSequence]; 375 updateBy2[iSequence]=0.0; 376 ClpSimplex::Status status = model_>getStatus(iSequence); 377 378 switch(status) { 379 380 case ClpSimplex::basic: 381 case ClpSimplex::isFixed: 382 break; 383 case ClpSimplex::isFree: 384 case ClpSimplex::superBasic: 385 if (fabs(value)>FREE_ACCEPT*tolerance) { 386 numberWanted; 387 // we are going to bias towards free (but only if reasonable) 388 if (fabs(value*FREE_BIAS)>bestDj) { 389 // check flagged variable and correct dj 390 if (!model_>flagged(iSequence)) { 391 bestDj=value*FREE_BIAS; 392 dj[iSequence]=value; 393 bestSequence = iSequence; 394 } else { 395 // just to make sure we don't exit before got something 396 numberWanted++; 397 } 398 } 399 } 400 break; 401 case ClpSimplex::atUpperBound: 402 if (value>tolerance) { 403 numberWanted; 404 if (value>bestDj) { 405 // check flagged variable and correct dj 406 if (!model_>flagged(iSequence)) { 407 bestDj=value; 408 dj[iSequence]=value; 409 bestSequence = iSequence; 410 } else { 411 // just to make sure we don't exit before got something 412 numberWanted++; 413 } 414 } 415 } 416 break; 417 case ClpSimplex::atLowerBound: 418 if (value<tolerance) { 419 numberWanted; 420 if (value>bestDj) { 421 // check flagged variable and correct dj 422 if (!model_>flagged(iSequence)) { 423 bestDj=value; 424 dj[iSequence]=value; 425 bestSequence = iSequence; 426 } else { 427 // just to make sure we don't exit before got something 428 numberWanted++; 429 } 430 } 431 } 432 break; 433 } 434 if (!numberWanted) { 435 // erase rest 436 for (;j<number;j++) { 437 int iSequence = index[j]; 438 updateBy2[iSequence]=0.0; 439 } 440 break; 441 } 442 } 443 if (!numberWanted) 444 break; 445 } 446 if (!numberWanted) 447 break; 448 } 449 spareColumn1>setNumElements(0); 450 spareColumn2>setNumElements(0); 451 updates>setDenseVector(saveDense); 452 /*if (model_>numberIterations()%100==0) 453 printf("%d best %g\n",bestSequence,bestDj);*/ 454 455 #ifdef CLP_DEBUG 456 if (bestSequence>=0) { 457 if (model_>getStatus(bestSequence)==ClpSimplex::atLowerBound) 458 assert(model_>reducedCost(bestSequence)<0.0); 459 if (model_>getStatus(bestSequence)==ClpSimplex::atUpperBound) 460 assert(model_>reducedCost(bestSequence)>0.0); 461 } 1876 thisWeight += pivotSquared * devex_ + pivot * modification; 1877 if (thisWeight<TRY_NORM) { 1878 if (mode_==1) { 1879 // steepest 1880 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 1881 } else { 1882 // exact 1883 thisWeight = referenceIn*pivotSquared; 1884 if (reference(iSequence)) 1885 thisWeight += 1.0; 1886 thisWeight = max(thisWeight,TRY_NORM); 1887 } 1888 } 1889 weight[iSequence] = thisWeight; 1890 } 1891 // restore outgoing weight 1892 if (sequenceOut>=0) 1893 weights_[sequenceOut]=outgoingWeight; 1894 alternateWeights_>clear(); 1895 spareColumn2>setNumElements(0); 1896 //#define SOME_DEBUG_1 1897 #ifdef SOME_DEBUG_1 1898 // check for accuracy 1899 int iCheck=229; 1900 //printf("weight for iCheck is %g\n",weights_[iCheck]); 1901 int numberRows = model_>numberRows(); 1902 //int numberColumns = model_>numberColumns(); 1903 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 1904 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 1905 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 1906 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 1907 } 462 1908 #endif 463 return bestSequence; 464 } else { 465 // Yes initialize 466 pivotSequence_=1; 467 pivotRow=1; 468 numberSwitched_++; 469 // Redo dualsolution 470 model_>computeDuals(NULL); 471 saveWeights(model_,4); 472 printf("switching %d nel ratio %g\n",numberElements,ratio); 473 updates>clear(); 474 } 475 } 476 } 1909 updates>setNumElements(0); 1910 spareColumn1>setNumElements(0); 1911 } 1912 // Returns pivot column, 1 if none 1913 int 1914 ClpPrimalColumnSteepest::pivotColumnOldMethod(CoinIndexedVector * updates, 1915 CoinIndexedVector * spareRow1, 1916 CoinIndexedVector * spareRow2, 1917 CoinIndexedVector * spareColumn1, 1918 CoinIndexedVector * spareColumn2) 1919 { 1920 assert(model_); 1921 int iSection,j; 1922 int number=0; 1923 int * index; 1924 double * updateBy; 1925 double * reducedCost; 1926 // dj could be very small (or even zero  take care) 1927 double dj = model_>dualIn(); 1928 double tolerance=model_>currentDualTolerance(); 1929 // we can't really trust infeasibilities if there is dual error 1930 // this coding has to mimic coding in checkDualSolution 1931 double error = min(1.0e3,model_>largestDualError()); 1932 // allow tolerance at least slightly bigger than standard 1933 tolerance = tolerance + error; 1934 int pivotRow = model_>pivotRow(); 477 1935 int anyUpdates; 478 1936 double * infeas = infeasible_>denseVector(); 479 1937 1938 // Local copy of mode so can decide what to do 1939 int switchType; 1940 if (mode_==4) 1941 switchType = 5numberSwitched_; 1942 else 1943 switchType = mode_; 1944 /* switchType  1945 0  all exact devex 1946 1  all steepest 1947 2  some exact devex 1948 3  auto some exact devex 1949 4  devex 1950 5  dantzig 1951 */ 480 1952 if (updates>getNumElements()) { 481 1953 // would have to have two goes for devex, three for steepest … … 854 2326 } 855 2327 } 2328 // See what sort of pricing 2329 int numberWanted=0; 2330 number = infeasible_>getNumElements(); 2331 int numberColumns = model_>numberColumns(); 2332 if (switchType==5) { 2333 // we can zero out 2334 updates>clear(); 2335 spareColumn1>clear(); 2336 pivotSequence_=1; 2337 pivotRow=1; 2338 // See if to switch 2339 int numberElements = model_>factorization()>numberElements(); 2340 int numberRows = model_>numberRows(); 2341 // ratio is done on number of columns here 2342 //double ratio = (double) numberElements/(double) numberColumns; 2343 double ratio = (double) numberElements/(double) numberRows; 2344 //double ratio = (double) numberElements/(double) model_>clpMatrix()>getNumElements(); 2345 if (ratio<0.1) { 2346 numberWanted = max(100,number/200); 2347 } else if (ratio<0.3) { 2348 numberWanted = max(500,number/40); 2349 } else if (ratio<0.5) { 2350 numberWanted = max(2000,number/10); 2351 numberWanted = max(numberWanted,numberColumns/30); 2352 } else { 2353 switchType=4; 2354 // initialize 2355 numberSwitched_++; 2356 // Make sure will redo 2357 delete [] weights_; 2358 weights_=NULL; 2359 saveWeights(model_,4); 2360 printf("switching to devex %d nel ratio %g\n",numberElements,ratio); 2361 updates>clear(); 2362 } 2363 if (model_>numberIterations()%1000==0) 2364 printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted); 2365 } 2366 if(switchType==4) { 2367 // Still in devex mode 2368 int numberElements = model_>factorization()>numberElements(); 2369 int numberRows = model_>numberRows(); 2370 // ratio is done on number of rows here 2371 double ratio = (double) numberElements/(double) numberRows; 2372 // Go to steepest if lot of iterations? 2373 if (ratio<1.0) { 2374 numberWanted = max(2000,number/20); 2375 } else if (ratio<5.0) { 2376 numberWanted = max(2000,number/10); 2377 numberWanted = max(numberWanted,numberColumns/20); 2378 } else { 2379 // we can zero out 2380 updates>clear(); 2381 spareColumn1>clear(); 2382 switchType=3; 2383 // initialize 2384 pivotSequence_=1; 2385 pivotRow=1; 2386 numberSwitched_++; 2387 // Make sure will redo 2388 delete [] weights_; 2389 weights_=NULL; 2390 saveWeights(model_,4); 2391 printf("switching to exact %d nel ratio %g\n",numberElements,ratio); 2392 updates>clear(); 2393 } 2394 if (model_>numberIterations()%1000==0) 2395 printf("numels %d ratio %g wanted %d\n",numberElements,ratio,numberWanted); 2396 } 2397 if (switchType<4) { 2398 if (switchType<2 ) { 2399 numberWanted = number+1; 2400 } else if (switchType==2) { 2401 numberWanted = max(2000,number/8); 2402 } else { 2403 int numberElements = model_>factorization()>numberElements(); 2404 double ratio = (double) numberElements/(double) model_>numberRows(); 2405 if (ratio<1.0) { 2406 numberWanted = max(2000,number/20); 2407 } else if (ratio<5.0) { 2408 numberWanted = max(2000,number/10); 2409 numberWanted = max(numberWanted,numberColumns/20); 2410 } else if (ratio<10.0) { 2411 numberWanted = max(2000,number/8); 2412 numberWanted = max(numberWanted,numberColumns/20); 2413 } else { 2414 ratio = number * (ratio/80.0); 2415 if (ratio>number) { 2416 numberWanted=number+1; 2417 } else { 2418 numberWanted = max(2000,(int) ratio); 2419 numberWanted = max(numberWanted,numberColumns/10); 2420 } 2421 } 2422 } 2423 } 856 2424 // for weights update we use pivotSequence 857 2425 pivotRow = pivotSequence_; … … 883 2451 updates,spareColumn2,spareColumn1); 884 2452 } 885 // now update weight update array 886 model_>factorization()>updateColumnTranspose(spareRow2, 887 alternateWeights_); 2453 double * weight; 888 2454 double * other = alternateWeights_>denseVector(); 889 double * weight;890 2455 int numberColumns = model_>numberColumns(); 891 2456 double scaleFactor = 1.0/dj; // as formula is with 1.0 … … 895 2460 updateBy = updates>denseVector(); 896 2461 weight = weights_+numberColumns; 897 898 for (j=0;j<number;j++) {899 int iSequence = index[j];900 double thisWeight = weight[iSequence];901 // row has 1902 double pivot = updateBy[iSequence]*scaleFactor;903 updateBy[iSequence]=0.0;904 double modification = other[iSequence];905 double pivotSquared = pivot * pivot;906 2462 907 thisWeight += pivotSquared * devex_ + pivot * modification; 908 if (thisWeight<TRY_NORM) { 909 if (switchType==1) { 910 // steepest 911 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 912 } else { 913 // exact 914 thisWeight = referenceIn*pivotSquared; 915 if (reference(iSequence+numberColumns)) 916 thisWeight += 1.0; 917 thisWeight = max(thisWeight,TRY_NORM); 918 } 919 } 920 weight[iSequence] = thisWeight; 2463 if (switchType<4) { 2464 // Exact 2465 // now update weight update array 2466 model_>factorization()>updateColumnTranspose(spareRow2, 2467 alternateWeights_); 2468 for (j=0;j<number;j++) { 2469 int iSequence = index[j]; 2470 double thisWeight = weight[iSequence]; 2471 // row has 1 2472 double pivot = updateBy[iSequence]*scaleFactor; 2473 updateBy[iSequence]=0.0; 2474 double modification = other[iSequence]; 2475 double pivotSquared = pivot * pivot; 2476 2477 thisWeight += pivotSquared * devex_ + pivot * modification; 2478 if (thisWeight<TRY_NORM) { 2479 if (switchType==1) { 2480 // steepest 2481 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 2482 } else { 2483 // exact 2484 thisWeight = referenceIn*pivotSquared; 2485 if (reference(iSequence+numberColumns)) 2486 thisWeight += 1.0; 2487 thisWeight = max(thisWeight,TRY_NORM); 2488 } 2489 } 2490 weight[iSequence] = thisWeight; 2491 } 2492 } else if (switchType==4) { 2493 // Devex 2494 assert (devex_>0.0); 2495 for (j=0;j<number;j++) { 2496 int iSequence = index[j]; 2497 double thisWeight = weight[iSequence]; 2498 // row has 1 2499 double pivot = updateBy[iSequence]*scaleFactor; 2500 updateBy[iSequence]=0.0; 2501 double value = pivot * pivot*devex_; 2502 if (reference(iSequence+numberColumns)) 2503 value += 1.0; 2504 weight[iSequence] = max(0.99*thisWeight,value); 2505 } 921 2506 } 922 2507 … … 929 2514 index = spareColumn1>getIndices(); 930 2515 updateBy = spareColumn1>denseVector(); 931 // get subset which have nonzero tableau elements 932 model_>clpMatrix()>subsetTransposeTimes(model_,alternateWeights_, 933 spareColumn1, 934 spareColumn2); 935 double * updateBy2 = spareColumn2>denseVector(); 936 for (j=0;j<number;j++) { 937 int iSequence = index[j]; 938 double thisWeight = weight[iSequence]; 939 double pivot = updateBy[iSequence]*scaleFactor; 940 updateBy[iSequence]=0.0; 941 double modification = updateBy2[iSequence]; 942 updateBy2[iSequence]=0.0; 943 double pivotSquared = pivot * pivot; 944 945 thisWeight += pivotSquared * devex_ + pivot * modification; 946 if (thisWeight<TRY_NORM) { 947 if (switchType==1) { 948 // steepest 949 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 950 } else { 951 // exact 952 thisWeight = referenceIn*pivotSquared; 953 if (reference(iSequence)) 954 thisWeight += 1.0; 955 thisWeight = max(thisWeight,TRY_NORM); 956 } 957 } 958 weight[iSequence] = thisWeight; 2516 if (switchType<4) { 2517 // Exact 2518 // get subset which have nonzero tableau elements 2519 model_>clpMatrix()>subsetTransposeTimes(model_,alternateWeights_, 2520 spareColumn1, 2521 spareColumn2); 2522 double * updateBy2 = spareColumn2>denseVector(); 2523 for (j=0;j<number;j++) { 2524 int iSequence = index[j]; 2525 double thisWeight = weight[iSequence]; 2526 double pivot = updateBy[iSequence]*scaleFactor; 2527 updateBy[iSequence]=0.0; 2528 double modification = updateBy2[iSequence]; 2529 updateBy2[iSequence]=0.0; 2530 double pivotSquared = pivot * pivot; 2531 2532 thisWeight += pivotSquared * devex_ + pivot * modification; 2533 if (thisWeight<TRY_NORM) { 2534 if (switchType==1) { 2535 // steepest 2536 thisWeight = max(TRY_NORM,ADD_ONE+pivotSquared); 2537 } else { 2538 // exact 2539 thisWeight = referenceIn*pivotSquared; 2540 if (reference(iSequence)) 2541 thisWeight += 1.0; 2542 thisWeight = max(thisWeight,TRY_NORM); 2543 } 2544 } 2545 weight[iSequence] = thisWeight; 2546 } 2547 } else if (switchType==4) { 2548 // Devex 2549 for (j=0;j<number;j++) { 2550 int iSequence = index[j]; 2551 double thisWeight = weight[iSequence]; 2552 // row has 1 2553 double pivot = updateBy[iSequence]*scaleFactor; 2554 updateBy[iSequence]=0.0; 2555 double value = pivot * pivot*devex_; 2556 if (reference(iSequence)) 2557 value += 1.0; 2558 weight[iSequence] = max(0.99*thisWeight,value); 2559 } 959 2560 } 960 2561 // restore outgoing weight … … 963 2564 alternateWeights_>clear(); 964 2565 spareColumn2>setNumElements(0); 2566 //#define SOME_DEBUG_1 965 2567 #ifdef SOME_DEBUG_1 966 2568 // check for accuracy 967 2569 int iCheck=229; 968 printf("weight for iCheck is %g\n",weights_[iCheck]); 969 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 970 !model_>getStatus(iCheck)!=isFixed) 971 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 2570 //printf("weight for iCheck is %g\n",weights_[iCheck]); 2571 int numberRows = model_>numberRows(); 2572 //int numberColumns = model_>numberColumns(); 2573 for (iCheck=0;iCheck<numberRows+numberColumns;iCheck++) { 2574 if (model_>getStatus(iCheck)!=ClpSimplex::basic&& 2575 !model_>getStatus(iCheck)!=ClpSimplex::isFixed) 2576 checkAccuracy(iCheck,1.0e1,updates,spareRow2); 2577 } 972 2578 #endif 973 2579 updates>setNumElements(0); … … 1008 2614 tolerance *= tolerance; // as we are using squares 1009 2615 1010 int numberWanted;1011 if (switchType<2 ) {1012 numberWanted = number+1;1013 } else if (switchType==2) {1014 numberWanted = max(2000,number/8);1015 } else {1016 int numberElements = model_>factorization()>numberElements();1017 double ratio = (double) numberElements/(double) model_>numberRows();1018 numberWanted = max(2000,number/8);1019 if (ratio<1.0) {1020 numberWanted = max(2000,number/20);1021 } else if (ratio>10.0) {1022 ratio = number * (ratio/80.0);1023 if (ratio>number)1024 numberWanted=number+1;1025 else1026 numberWanted = max(2000,(int) ratio);1027 }1028 }1029 2616 int iPass; 1030 2617 // Setup two passes … … 1039 2626 for (iPass=0;iPass<2;iPass++) { 1040 2627 int end = start[2*iPass+1]; 1041 for (i=start[2*iPass];i<end;i++) { 1042 iSequence = index[i]; 1043 double value = infeas[iSequence]; 1044 if (value>tolerance) { 1045 double weight = weights_[iSequence]; 1046 //weight=1.0; 1047 if (value>bestDj*weight) { 1048 // check flagged variable and correct dj 1049 if (!model_>flagged(iSequence)) { 1050 bestDj=value/weight; 1051 bestSequence = iSequence; 1052 } else { 1053 // just to make sure we don't exit before got something 1054 numberWanted++; 2628 if (switchType<5) { 2629 for (i=start[2*iPass];i<end;i++) { 2630 iSequence = index[i]; 2631 double value = infeas[iSequence]; 2632 if (value>tolerance) { 2633 double weight = weights_[iSequence]; 2634 //weight=1.0; 2635 if (value>bestDj*weight) { 2636 // check flagged variable and correct dj 2637 if (!model_>flagged(iSequence)) { 2638 bestDj=value/weight; 2639 bestSequence = iSequence; 2640 } else { 2641 // just to make sure we don't exit before got something 2642 numberWanted++; 2643 } 1055 2644 } 1056 2645 } 1057 } 1058 numberWanted; 1059 if (!numberWanted) 1060 break; 2646 numberWanted; 2647 if (!numberWanted) 2648 break; 2649 } 2650 } else { 2651 // Dantzig 2652 for (i=start[2*iPass];i<end;i++) { 2653 iSequence = index[i]; 2654 double value = infeas[iSequence]; 2655 if (value>tolerance) { 2656 if (value>bestDj) { 2657 // check flagged variable and correct dj 2658 if (!model_>flagged(iSequence)) { 2659 bestDj=value; 2660 bestSequence = iSequence; 2661 } else { 2662 // just to make sure we don't exit before got something 2663 numberWanted++; 2664 } 2665 } 2666 } 2667 numberWanted; 2668 if (!numberWanted) 2669 break; 2670 } 1061 2671 } 1062 2672 if (!numberWanted) … … 1092 2702 if (mode==1&&!weights_) 1093 2703 numberSwitched_=0; // Reset 1094 if(!numberSwitched_)1095 return;1096 2704 } 1097 2705 // alternateWeights_ is defined as indexed but is treated oddly … … 1341 2949 alternateWeights_>setNumElements(number); 1342 2950 } else { 1343 for (i=0;i<number;i++) { 1344 int iRow = which[i]; 1345 int iPivot = pivotVariable[iRow]; 1346 if (reference(iPivot)) { 1347 devex_ += work[iRow]*work[iRow]; 1348 newWork[iRow] = 2.0*work[iRow]; 1349 newWhich[newNumber++]=iRow; 1350 } 1351 } 1352 if (!newWork[pivotRow]) 1353 newWhich[newNumber++]=pivotRow; // add if not already in 1354 newWork[pivotRow] = 2.0*max(devex_,0.0); 2951 if (mode_!=4numberSwitched_>1) { 2952 for (i=0;i<number;i++) { 2953 int iRow = which[i]; 2954 int iPivot = pivotVariable[iRow]; 2955 if (reference(iPivot)) { 2956 devex_ += work[iRow]*work[iRow]; 2957 newWork[iRow] = 2.0*work[iRow]; 2958 newWhich[newNumber++]=iRow; 2959 } 2960 } 2961 if (!newWork[pivotRow]) 2962 newWhich[newNumber++]=pivotRow; // add if not already in 2963 newWork[pivotRow] = 2.0*max(devex_,0.0); 2964 } else { 2965 for (i=0;i<number;i++) { 2966 int iRow = which[i]; 2967 int iPivot = pivotVariable[iRow]; 2968 if (reference(iPivot)) 2969 devex_ += work[iRow]*work[iRow]; 2970 } 2971 } 1355 2972 if (reference(sequenceIn)) { 1356 2973 devex_+=1.0; … … 1388 3005 printf("old weight %g, new %g\n",oldDevex,devex_); 1389 3006 #endif 1390 double check = max(devex_,oldDevex) ;;3007 double check = max(devex_,oldDevex)+0.1; 1391 3008 weights_[sequenceIn] = devex_; 1392 if ( fabs ( devex_  oldDevex ) > 1.0e1 * check ) { 3009 double testValue=0.1; 3010 if (mode_==4&&numberSwitched_==1) 3011 testValue = 0.5; 3012 if ( fabs ( devex_  oldDevex ) > testValue * check ) { 1393 3013 #ifdef CLP_DEBUG 1394 3014 if ((model_>messageHandler()>logLevel()&48)==16) 1395 3015 printf("old weight %g, new %g\n",oldDevex,devex_); 1396 3016 #endif 1397 if ( fabs ( devex_  oldDevex ) > 1.0e5 * check ) { 3017 //printf("old weight %g, new %g\n",oldDevex,devex_); 3018 testValue=0.99; 3019 if (mode_==1) 3020 testValue=1.01e1; // make unlikely to do if steepest 3021 else if (mode_==4&&numberSwitched_==1) 3022 testValue=0.9; 3023 double difference = fabs(devex_oldDevex); 3024 if ( difference > testValue * check ) { 1398 3025 // need to redo 1399 3026 model_>messageHandler()>message(CLP_INITIALIZE_STEEP,
Note: See TracChangeset
for help on using the changeset viewer.