Changeset 225 for trunk/ClpNetworkMatrix.cpp
 Timestamp:
 Oct 15, 2003 5:34:57 PM (16 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/ClpNetworkMatrix.cpp
r124 r225 15 15 #include "ClpPlusMinusOneMatrix.hpp" 16 16 #include "ClpMessage.hpp" 17 #include <iostream> 17 18 18 19 //############################################################################# … … 369 370 ClpPlusMinusOneMatrix* rowCopy = 370 371 dynamic_cast< ClpPlusMinusOneMatrix*>(model>rowCopy()); 371 if (numberInRowArray>0.3*numberRows!rowCopy) { 372 bool packed = rowArray>packedMode(); 373 double factor = 0.3; 374 // We may not want to do by row if there may be cache problems 375 int numberColumns = model>numberColumns(); 376 // It would be nice to find L2 cache size  for moment 512K 377 // Be slightly optimistic 378 if (numberColumns*sizeof(double)>1000000) { 379 if (numberRows*10<numberColumns) 380 factor=0.1; 381 else if (numberRows*4<numberColumns) 382 factor=0.15; 383 else if (numberRows*2<numberColumns) 384 factor=0.2; 385 //if (model>numberIterations()%50==0) 386 //printf("%d nonzero\n",numberInRowArray); 387 } 388 if (numberInRowArray>factor*numberRows!rowCopy) { 372 389 // do by column 373 390 int iColumn; 374 double * markVector = y>denseVector(); // probably empty391 assert (!y>getNumElements()); 375 392 CoinBigIndex j=0; 376 if (trueNetwork_) { 377 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 378 double value = markVector[iColumn]; 379 markVector[iColumn]=0.0; 380 int iRowM = indices_[j]; 381 int iRowP = indices_[j+1]; 382 value = scalar*pi[iRowM]; 383 value += scalar*pi[iRowP]; 384 if (fabs(value)>zeroTolerance) { 385 index[numberNonZero++]=iColumn; 386 array[iColumn]=value; 387 } 393 if (packed) { 394 // need to expand pi into y 395 assert(y>capacity()>=numberRows); 396 double * piOld = pi; 397 pi = y>denseVector(); 398 const int * whichRow = rowArray>getIndices(); 399 int i; 400 // modify pi so can collapse to one loop 401 for (i=0;i<numberInRowArray;i++) { 402 int iRow = whichRow[i]; 403 pi[iRow]=scalar*piOld[i]; 404 } 405 if (trueNetwork_) { 406 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 407 double value = 0.0; 408 int iRowM = indices_[j]; 409 int iRowP = indices_[j+1]; 410 value = pi[iRowM]; 411 value += pi[iRowP]; 412 if (fabs(value)>zeroTolerance) { 413 array[numberNonZero]=value; 414 index[numberNonZero++]=iColumn; 415 } 416 } 417 } else { 418 // skip negative rows 419 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 420 double value = 0.0; 421 int iRowM = indices_[j]; 422 int iRowP = indices_[j+1]; 423 if (iRowM>=0) 424 value = pi[iRowM]; 425 if (iRowP>=0) 426 value += pi[iRowP]; 427 if (fabs(value)>zeroTolerance) { 428 array[numberNonZero]=value; 429 index[numberNonZero++]=iColumn; 430 } 431 } 432 } 433 for (i=0;i<numberInRowArray;i++) { 434 int iRow = whichRow[i]; 435 pi[iRow]=0.0; 388 436 } 389 437 } else { 390 // skip negative rows 391 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 392 double value = markVector[iColumn]; 393 markVector[iColumn]=0.0; 394 int iRowM = indices_[j]; 395 int iRowP = indices_[j+1]; 396 if (iRowM>=0) 438 if (trueNetwork_) { 439 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 440 double value = 0.0; 441 int iRowM = indices_[j]; 442 int iRowP = indices_[j+1]; 397 443 value = scalar*pi[iRowM]; 398 if (iRowP>=0)399 444 value += scalar*pi[iRowP]; 400 if (fabs(value)>zeroTolerance) { 401 index[numberNonZero++]=iColumn; 402 array[iColumn]=value; 445 if (fabs(value)>zeroTolerance) { 446 index[numberNonZero++]=iColumn; 447 array[iColumn]=value; 448 } 449 } 450 } else { 451 // skip negative rows 452 for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) { 453 double value = 0.0; 454 int iRowM = indices_[j]; 455 int iRowP = indices_[j+1]; 456 if (iRowM>=0) 457 value = scalar*pi[iRowM]; 458 if (iRowP>=0) 459 value += scalar*pi[iRowP]; 460 if (fabs(value)>zeroTolerance) { 461 index[numberNonZero++]=iColumn; 462 array[iColumn]=value; 463 } 403 464 } 404 465 } 405 466 } 406 467 columnArray>setNumElements(numberNonZero); 407 y>setNumElements(0);408 468 } else { 409 469 // do by row … … 430 490 int numberToDo = y>getNumElements(); 431 491 const int * which = y>getIndices(); 432 if (trueNetwork_) { 433 for (jColumn=0;jColumn<numberToDo;jColumn++) { 434 int iColumn = which[jColumn]; 435 double value = 0.0; 436 CoinBigIndex j=iColumn<<1; 437 int iRowM = indices_[j]; 438 int iRowP = indices_[j+1]; 439 value = pi[iRowM]; 440 value += pi[iRowP]; 441 if (fabs(value)>zeroTolerance) { 442 index[numberNonZero++]=iColumn; 443 array[iColumn]=value; 444 } 492 bool packed = rowArray>packedMode(); 493 if (packed) { 494 // Must line up with y 495 // need to expand pi into y 496 int numberInRowArray = rowArray>getNumElements(); 497 assert(y>capacity()>=model>numberRows()); 498 double * piOld = pi; 499 pi = y>denseVector(); 500 const int * whichRow = rowArray>getIndices(); 501 int i; 502 for (i=0;i<numberInRowArray;i++) { 503 int iRow = whichRow[i]; 504 pi[iRow]=piOld[i]; 505 } 506 if (trueNetwork_) { 507 for (jColumn=0;jColumn<numberToDo;jColumn++) { 508 int iColumn = which[jColumn]; 509 double value = 0.0; 510 CoinBigIndex j=iColumn<<1; 511 int iRowM = indices_[j]; 512 int iRowP = indices_[j+1]; 513 value = pi[iRowM]; 514 value += pi[iRowP]; 515 array[jColumn]=value; 516 } 517 } else { 518 // skip negative rows 519 for (jColumn=0;jColumn<numberToDo;jColumn++) { 520 int iColumn = which[jColumn]; 521 double value = 0.0; 522 CoinBigIndex j=iColumn<<1; 523 int iRowM = indices_[j]; 524 int iRowP = indices_[j+1]; 525 if (iRowM>=0) 526 value = pi[iRowM]; 527 if (iRowP>=0) 528 value += pi[iRowP]; 529 array[jColumn]=value; 530 } 531 } 532 for (i=0;i<numberInRowArray;i++) { 533 int iRow = whichRow[i]; 534 pi[iRow]=0.0; 445 535 } 446 536 } else { 447 // skip negative rows 448 for (jColumn=0;jColumn<numberToDo;jColumn++) { 449 int iColumn = which[jColumn]; 450 double value = 0.0; 451 CoinBigIndex j=iColumn<<1; 452 int iRowM = indices_[j]; 453 int iRowP = indices_[j+1]; 454 if (iRowM>=0) 537 if (trueNetwork_) { 538 for (jColumn=0;jColumn<numberToDo;jColumn++) { 539 int iColumn = which[jColumn]; 540 double value = 0.0; 541 CoinBigIndex j=iColumn<<1; 542 int iRowM = indices_[j]; 543 int iRowP = indices_[j+1]; 455 544 value = pi[iRowM]; 456 if (iRowP>=0)457 545 value += pi[iRowP]; 458 if (fabs(value)>zeroTolerance) { 459 index[numberNonZero++]=iColumn; 460 array[iColumn]=value; 546 if (fabs(value)>zeroTolerance) { 547 index[numberNonZero++]=iColumn; 548 array[iColumn]=value; 549 } 550 } 551 } else { 552 // skip negative rows 553 for (jColumn=0;jColumn<numberToDo;jColumn++) { 554 int iColumn = which[jColumn]; 555 double value = 0.0; 556 CoinBigIndex j=iColumn<<1; 557 int iRowM = indices_[j]; 558 int iRowP = indices_[j+1]; 559 if (iRowM>=0) 560 value = pi[iRowM]; 561 if (iRowP>=0) 562 value += pi[iRowP]; 563 if (fabs(value)>zeroTolerance) { 564 index[numberNonZero++]=iColumn; 565 array[iColumn]=value; 566 } 461 567 } 462 568 } … … 542 648 return numberElements; 543 649 } 650 /* If element NULL returns number of elements in column part of basis, 651 If not NULL fills in as well */ 652 CoinBigIndex 653 ClpNetworkMatrix::fillBasis(const ClpSimplex * model, 654 const int * whichColumn, 655 int numberBasic, 656 int numberColumnBasic, 657 int * indexRowU, int * indexColumnU, 658 double * elementU) const 659 { 660 int i; 661 CoinBigIndex numberElements=0; 662 if (elementU!=NULL) { 663 if (trueNetwork_) { 664 for (i=0;i<numberColumnBasic;i++) { 665 int iColumn = whichColumn[i]; 666 CoinBigIndex j=iColumn<<1; 667 int iRowM = indices_[j]; 668 int iRowP = indices_[j+1]; 669 indexRowU[numberElements]=iRowM; 670 indexColumnU[numberElements]=numberBasic; 671 elementU[numberElements]=1.0; 672 indexRowU[numberElements+1]=iRowP; 673 indexColumnU[numberElements+1]=numberBasic; 674 elementU[numberElements+1]=1.0; 675 numberElements+=2; 676 numberBasic++; 677 } 678 } else { 679 for (i=0;i<numberColumnBasic;i++) { 680 int iColumn = whichColumn[i]; 681 CoinBigIndex j=iColumn<<1; 682 int iRowM = indices_[j]; 683 int iRowP = indices_[j+1]; 684 if (iRowM>=0) { 685 indexRowU[numberElements]=iRowM; 686 indexColumnU[numberElements]=numberBasic; 687 elementU[numberElements++]=1.0; 688 } 689 if (iRowP>=0) { 690 indexRowU[numberElements]=iRowP; 691 indexColumnU[numberElements]=numberBasic; 692 elementU[numberElements++]=1.0; 693 } 694 numberBasic++; 695 } 696 } 697 } else { 698 if (trueNetwork_) { 699 numberElements = 2*numberColumnBasic; 700 } else { 701 for (i=0;i<numberColumnBasic;i++) { 702 int iColumn = whichColumn[i]; 703 CoinBigIndex j=iColumn<<1; 704 int iRowM = indices_[j]; 705 int iRowP = indices_[j+1]; 706 if (iRowM>=0) 707 numberElements ++; 708 if (iRowP>=0) 709 numberElements ++; 710 } 711 } 712 } 713 return numberElements; 714 } 544 715 /* Unpacks a column into an CoinIndexedvector 545 Note that model is NOT const. Bounds and objective could 546 be modified if doing column generation */ 716 */ 547 717 void 548 718 ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray, … … 557 727 rowArray>add(iRowP,1.0); 558 728 } 729 /* Unpacks a column into an CoinIndexedvector 730 ** in packed foramt 731 Note that model is NOT const. Bounds and objective could 732 be modified if doing column generation (just for this variable) */ 733 void 734 ClpNetworkMatrix::unpackPacked(ClpSimplex * model, 735 CoinIndexedVector * rowArray, 736 int iColumn) const 737 { 738 int * index = rowArray>getIndices(); 739 double * array = rowArray>denseVector(); 740 int number = 0; 741 CoinBigIndex j=iColumn<<1; 742 int iRowM = indices_[j]; 743 int iRowP = indices_[j+1]; 744 if (iRowM>=0) { 745 array[number]=1.0; 746 index[number++]=iRowM; 747 } 748 if (iRowP>=0) { 749 array[number]=1.0; 750 index[number++]=iRowP; 751 } 752 rowArray>setNumElements(number); 753 rowArray>setPackedMode(true); 754 } 559 755 /* Adds multiple of a column into an CoinIndexedvector 560 756 You can use quickAdd to add to vector */ … … 631 827 void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 632 828 { 829 std::cerr<<"deleteCols not implemented in ClpNetworkMatrix"<<std::endl; 633 830 abort(); 634 831 } … … 636 833 void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 637 834 { 835 std::cerr<<"deleteRows not implemented in ClpNetworkMatrix"<<std::endl; 638 836 abort(); 639 837 } 838 /* Given positive integer weights for each row fills in sum of weights 839 for each column (and slack). 840 Returns weights vector 841 */ 842 CoinBigIndex * 843 ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const 844 { 845 int numberRows = model>numberRows(); 846 int numberColumns =model>numberColumns(); 847 int number = numberRows+numberColumns; 848 CoinBigIndex * weights = new CoinBigIndex[number]; 849 int i; 850 for (i=0;i<numberColumns;i++) { 851 CoinBigIndex j=i<<1; 852 CoinBigIndex count=0; 853 int iRowM = indices_[j]; 854 int iRowP = indices_[j+1]; 855 if (iRowM>=0) { 856 count += inputWeights[iRowM]; 857 } 858 if (iRowP>=0) { 859 count += inputWeights[iRowP]; 860 } 861 weights[i]=count; 862 } 863 for (i=0;i<numberRows;i++) { 864 weights[i+numberColumns]=inputWeights[i]; 865 } 866 return weights; 867 } 868 /* Returns largest and smallest elements of both signs. 869 Largest refers to largest absolute value. 870 */ 871 void 872 ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative, 873 double & smallestPositive, double & largestPositive) 874 { 875 smallestNegative=1.0; 876 largestNegative=1.0; 877 smallestPositive=1.0; 878 largestPositive=1.0; 879 }
Note: See TracChangeset
for help on using the changeset viewer.