Changeset 364
 Timestamp:
 May 11, 2004 3:36:06 PM (15 years ago)
 Location:
 trunk
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

trunk/ClpPredictorCorrector.cpp
r362 r364 99 99 double * savePrimal=NULL; 100 100 int numberTotal = numberRows_+numberColumns_; 101 // Extra regions for centering 102 double * saveWeights = new double[numberTotal]; 103 double * saveZ = new double[numberTotal]; 104 double * saveW = new double[numberTotal]; 101 105 while (problemStatus_<0) { 102 106 if (numberIterations_==xxxxxx) { … … 309 313 bool goodMove=false; 310 314 double bestNextGap=COIN_DBL_MAX; 315 double nextGap=bestNextGap; 316 int nextNumber=0; 311 317 worstDirectionAccuracy_=0.0; 312 318 while (!goodMove) { … … 370 376 worstDirectionAccuracy_=directionAccuracy; 371 377 } 372 findStepLength(phase); 373 int nextNumber; //number of complementarity pairs 374 double nextGap=complementarityGap(nextNumber,1); 378 findStepLength(phase,NULL); 379 nextGap=complementarityGap(nextNumber,1); 375 380 376 381 bestNextGap=max(nextGap,0.8*complementarityGap_); … … 473 478 } 474 479 if (goodMove) { 475 findStepLength(phase );480 findStepLength(phase,NULL); 476 481 // just for debug 477 int nextNumber; 478 complementarityGap(nextNumber,1); 482 nextGap = complementarityGap(nextNumber,1); 479 483 if (numberIterations_>=77) { 480 484 goodMove=checkGoodMove(true,bestNextGap); … … 507 511 testValue=maximumRHSError_; 508 512 } 509 findStepLength(2 );513 findStepLength(2,NULL); 510 514 // just for debug 511 int nextNumber; 512 complementarityGap(nextNumber,2); 515 nextGap=complementarityGap(nextNumber,2); 513 516 } 514 517 } /* endwhile */ … … 554 557 } 555 558 } 559 #if 0 560 // See if we can do centering steps 561 memcpy(saveWeights,weights_,numberTotal*sizeof(double)); 562 memcpy(saveZ,deltaZ_,numberTotal*sizeof(double)); 563 memcpy(saveW,deltaW_,numberTotal*sizeof(double)); 564 double savePrimalStep = actualPrimalStep_; 565 double saveDualStep = actualDualStep_; 566 double saveMu = mu_; 567 mu_=nextGap / (double) nextNumber; 568 setupForSolve(3); 569 findDirectionVector(3); 570 memcpy(deltaZ_,saveZ,numberTotal*sizeof(double)); 571 memcpy(deltaW_,saveW,numberTotal*sizeof(double)); 572 findStepLength(3,saveWeights); 573 double xGap = complementarityGap(nextNumber,3); 574 #if 0 575 goodMove=checkGoodMove(true,bestNextGap); 576 #endif 577 if (xGap>nextGap) { 578 mu_=saveMu; 579 actualPrimalStep_ = savePrimalStep; 580 actualDualStep_ = saveDualStep; 581 memcpy(weights_,saveWeights,numberTotal*sizeof(double)); 582 memcpy(deltaZ_,saveZ,numberTotal*sizeof(double)); 583 memcpy(deltaW_,saveW,numberTotal*sizeof(double)); 584 } 585 #endif 556 586 numberFixed=updateSolution(); 557 587 numberFixedTotal+=numberFixed; 558 588 } /* endwhile */ 589 delete [] saveWeights; 590 delete [] saveZ; 591 delete [] saveW; 559 592 if (savePi) { 560 593 //std::cout<<"Restoring from iteration "<<saveIteration<<std::endl; … … 595 628 // 1 corrector 596 629 // 2 primal dual 597 double ClpPredictorCorrector::findStepLength(const int phase )630 double ClpPredictorCorrector::findStepLength(const int phase,const double * oldWeight) 598 631 { 599 632 double directionNorm=0.0; … … 870 903 } 871 904 break; 872 } 873 printf("step  phase %d, norm %g, dual step %g, primal step %g\n", 874 phase,directionNorm,maximumPrimalStep,maximumDualStep); 905 case 3: 906 //centering step 907 for (int iColumn=0;iColumn<numberTotal;iColumn++) { 908 if (!flagged(iColumn)) { 909 // old deltas from previous step 910 double z1=work1[iColumn]; 911 double w1=work2[iColumn]; 912 double work3Value=0.0; 913 double work4Value=0.0; 914 double directionElement=weights[iColumn]; 915 double oldElement = oldWeight[iColumn]; 916 double newElement = oldElement+directionElement; 917 weights[iColumn]=newElement; 918 double value=primal[iColumn]; 919 if (directionNorm<fabs(directionElement)) { 920 directionNorm=fabs(directionElement); 921 } 922 if (lowerBound(iColumn) 923 upperBound(iColumn)) { 924 if (lowerBound(iColumn)) { 925 double gap=lowerSlack[iColumn]+extra; 926 double delta; 927 double zDelta; 928 if (!fakeLower(iColumn)) { 929 zDelta = lower[iColumn]+lowerSlack[iColumn] 930  value  directionElement; 931 } else { 932 zDelta =  directionElement; 933 } 934 delta = zDelta  oldElement; 935 z1+=(work3[iColumn]zVec[iColumn]*zDelta)/gap; 936 work3Value=delta; 937 if (lowerSlack[iColumn]<maximumPrimalStep*delta) { 938 maximumPrimalStep=lowerSlack[iColumn]/delta; 939 chosenPrimalSequence=iColumn; 940 } 941 if (zVec[iColumn]>tolerance) { 942 if (zVec[iColumn]<z1*maximumDualStep) { 943 maximumDualStep=zVec[iColumn]/z1; 944 chosenDualSequence=iColumn; 945 } 946 } 947 } 948 if (upperBound(iColumn)) { 949 double gap=upperSlack[iColumn]+extra; 950 double delta; 951 double zDelta; 952 if (!fakeUpper(iColumn)) { 953 zDelta = upper[iColumn]+upperSlack[iColumn] 954 + value + directionElement; 955 } else { 956 zDelta = directionElement; 957 } 958 delta = zDelta + oldElement; 959 //double delta = upper[iColumn]+upperSlack[iColumn] 960 //+ value + directionElement; 961 w1+=(work4[iColumn]+wVec[iColumn]*zDelta)/gap; 962 work4Value=delta; 963 if (upperSlack[iColumn]<maximumPrimalStep*delta) { 964 maximumPrimalStep=upperSlack[iColumn]/delta; 965 chosenPrimalSequence=iColumn; 966 } 967 if (wVec[iColumn]>tolerance) { 968 if (wVec[iColumn]<w1*maximumDualStep) { 969 maximumDualStep=wVec[iColumn]/w1; 970 chosenDualSequence=iColumn; 971 } 972 } 973 } 974 } else { 975 //free 976 double gap=fabs(value); 977 double multiplier=1.0/gap; 978 if (gap<1.0) { 979 multiplier=1,0; 980 } 981 z1+=multiplier*(work3[iColumn]directionElement*zVec[iColumn]); 982 w1+=multiplier*(work4[iColumn]+directionElement*wVec[iColumn]); 983 } 984 work1[iColumn]=z1; 985 work2[iColumn]=w1; 986 work3[iColumn]=work3Value; 987 work4[iColumn]=work4Value; 988 } else { 989 work1[iColumn]=0.0; 990 work2[iColumn]=0.0; 991 work3[iColumn]=0.0; 992 work4[iColumn]=0.0; 993 } 994 } 995 break; 996 } 997 //printf("step  phase %d, norm %g, dual step %g, primal step %g\n", 998 // phase,directionNorm,maximumPrimalStep,maximumDualStep); 875 999 actualPrimalStep_=stepLength_*maximumPrimalStep; 876 1000 if (phase>=0&&actualPrimalStep_>1.0) { … … 956 1080 double ClpPredictorCorrector::findDirectionVector(const int phase) 957 1081 { 1082 if (phase==3) { 1083 // centering step 1084 double * work2 = deltaW_; 1085 int numberTotal = numberRows_+numberColumns_; 1086 double * tempRegion = new double [numberRows_]; 1087 double * newError = new double [numberRows_]; 1088 // For KKT separate out 1089 #ifndef KKT 1090 int iColumn; 1091 for (iColumn=0;iColumn<numberTotal;iColumn++) 1092 weights_[iColumn] = work2[iColumn]; 1093 multiplyAdd(weights_+numberColumns_,numberRows_,1.0,tempRegion,0.0); 1094 matrix_>times(1.0,weights_,tempRegion); 1095 #else 1096 // regions in will be work2 and newError 1097 // regions out weights_ and tempRegion 1098 memset(newError,0,numberRows_*sizeof(double)); 1099 double * region1Save=NULL;//for refinement 1100 #endif 1101 #ifndef KKT 1102 double maximumRHS = maximumAbsElement(tempRegion,numberRows_); 1103 double scale=1.0; 1104 double unscale=1.0; 1105 if (maximumRHS>1.0e30) { 1106 if (maximumRHS<=0.5) { 1107 double factor=2.0; 1108 while (maximumRHS<=0.5) { 1109 maximumRHS*=factor; 1110 scale*=factor; 1111 } /* endwhile */ 1112 } else if (maximumRHS>=2.0&&maximumRHS<=COIN_DBL_MAX) { 1113 double factor=0.5; 1114 while (maximumRHS>=2.0) { 1115 maximumRHS*=factor; 1116 scale*=factor; 1117 } /* endwhile */ 1118 } 1119 unscale=diagonalScaleFactor_/scale; 1120 } else { 1121 //effectively zero 1122 scale=0.0; 1123 unscale=0.0; 1124 } 1125 //printf("putting scales to 1.0\n"); 1126 //scale=1.0; 1127 //unscale=1.0; 1128 multiplyAdd(NULL,numberRows_,0.0,tempRegion,scale); 1129 cholesky_>solve(tempRegion); 1130 multiplyAdd(NULL,numberRows_,0.0,tempRegion,unscale); 1131 //CoinZeroN(newError,numberRows_); 1132 multiplyAdd(tempRegion,numberRows_,1.0,weights_+numberColumns_,0.0); 1133 CoinZeroN(weights_,numberColumns_); 1134 matrix_>transposeTimes(1.0,tempRegion,weights_); 1135 //if flagged then entries zero so can do 1136 for (iColumn=0;iColumn<numberTotal;iColumn++) 1137 weights_[iColumn] = weights_[iColumn]*diagonal_[iColumn] 1138 work2[iColumn]; 1139 #else 1140 solveSystem(weights_, tempRegion, 1141 work2,newError,region1Save,regionSave,lastError>1.0e5); 1142 #endif 1143 multiplyAdd(weights_+numberColumns_,numberRows_,1.0,newError,0.0); 1144 matrix_>times(1.0,weights_,newError); 1145 delete [] newError; 1146 delete [] tempRegion; 1147 return 0.0; 1148 } 958 1149 double projectionTolerance=projectionTolerance_; 959 1150 //temporary … … 1622 1813 numberComplementarityPairs=1; 1623 1814 } 1624 printf("gap %g  smallest %g, largest %g, pairs %d\n",1625 1815 //printf("gap %g  smallest %g, largest %g, pairs %d\n", 1816 // gap,smallestGap,largestGap,numberComplementarityPairs); 1626 1817 return gap; 1627 1818 } … … 1644 1835 1645 1836 int iColumn; 1646 printf("phase %d in setupForSolve, mu %g\n",phase,mu_);1837 //printf("phase %d in setupForSolve, mu %g\n",phase,mu_); 1647 1838 switch (phase) { 1648 1839 case 0: … … 1739 1930 break; 1740 1931 case 3: 1741 for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) { 1742 if (!flagged(iColumn)) { 1743 double value= 0.0; 1744 if (lowerBound(iColumn)) { 1745 value=work3[iColumn]/(lowerSlack[iColumn]+extra); 1746 } 1747 if (upperBound(iColumn)) { 1748 value+=work4[iColumn]/(upperSlack[iColumn]+extra); 1749 } 1932 { 1933 double minBeta = 0.1*mu_; 1934 double maxBeta = 10.0*mu_; 1935 double * work1 = deltaZ_; 1936 double * weights = weights_; 1937 for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) { 1938 work3[iColumn]=0.0; 1939 work4[iColumn]=0.0; 1940 if (!flagged(iColumn)) { 1941 double value= 0.0; 1942 if (lowerBound(iColumn)) { 1943 double change; 1944 if (!fakeLower(iColumn)) { 1945 change =primal[iColumn]+weights[iColumn]lowerSlack[iColumn]lower[iColumn]; 1946 } else { 1947 change =weights[iColumn]; 1948 } 1949 double dualValue=zVec[iColumn]+actualDualStep_*work1[iColumn]; 1950 double primalValue=lowerSlack[iColumn]+actualPrimalStep_*change; 1951 double gapProduct=dualValue*primalValue; 1952 if (gapProduct<minBeta) { 1953 value= (minBetagapProduct)/(lowerSlack[iColumn]+extra); 1954 } else if (gapProduct>maxBeta) { 1955 value= max(maxBetagapProduct,maxBeta)/(lowerSlack[iColumn]+extra); 1956 } 1957 value=value; 1958 work3[iColumn]=value; 1959 } 1960 if (upperBound(iColumn)) { 1961 double change; 1962 if (!fakeUpper(iColumn)) { 1963 change =upper[iColumn]primal[iColumn]weights[iColumn]upperSlack[iColumn]; 1964 } else { 1965 change =weights[iColumn]; 1966 } 1967 double dualValue=wVec[iColumn]+actualDualStep_*work2[iColumn]; 1968 double primalValue=upperSlack[iColumn]+actualPrimalStep_*change; 1969 double gapProduct=dualValue*primalValue; 1970 if (gapProduct<minBeta) { 1971 value+= (minBetagapProduct)/(upperSlack[iColumn]+extra); 1972 } else if (gapProduct>maxBeta) { 1973 value+= max(maxBetagapProduct,maxBeta)/(upperSlack[iColumn]+extra); 1974 } 1975 value=value; 1976 work4[iColumn]=valuework3[iColumn]; 1977 } 1750 1978 #ifndef KKT 1751 work2[iColumn]=diagonal_[iColumn]*(dual[iColumn]+value);1979 work2[iColumn]=diagonal_[iColumn]*value; 1752 1980 #else 1753 work2[iColumn]=dual[iColumn]+value;1981 work2[iColumn]=value; 1754 1982 #endif 1755 } else { 1756 work2[iColumn]=0.0; 1757 } 1758 } 1983 } else { 1984 work2[iColumn]=0.0; 1985 } 1986 } 1987 } 1759 1988 break; 1760 1989 } /* endswitch */ 
trunk/include/ClpPredictorCorrector.hpp
r328 r364 44 44 // 1 corrector 45 45 // 2 primal dual 46 double findStepLength(const int phase );46 double findStepLength(const int phase,const double * oldWeight); 47 47 /// findDirectionVector. 48 48 double findDirectionVector(const int phase);
Note: See TracChangeset
for help on using the changeset viewer.