Changeset 1761
 Timestamp:
 Jul 6, 2011 12:06:24 PM (8 years ago)
 Location:
 trunk/Clp/src
 Files:

 9 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpEventHandler.hpp
r1665 r1761 34 34 enum Event { 35 35 endOfIteration = 100, // used to set secondary status 36 endOfFactorization, 36 endOfFactorization, // after gutsOfSolution etc 37 37 endOfValuesPass, 38 38 node, // for Cbc … … 47 47 presolveAfterFirstSolve, // ClpSolve presolve after solve 48 48 presolveAfterSolve, // ClpSolve presolve after solve 49 presolveEnd // ClpSolve presolve end 49 presolveEnd, // ClpSolve presolve end 50 goodFactorization, // before gutsOfSolution 51 complicatedPivotIn, // in modifyCoefficients 52 noCandidateInPrimal, // tentative end 53 looksEndInPrimal, // About to declare victory (or defeat) 54 endInPrimal, // Victory (or defeat) 55 complicatedPivotOut, // in modifyCoefficients 56 noCandidateInDual, // tentative end 57 looksEndInDual, // About to declare victory (or defeat) 58 endInDual, // Victory (or defeat) 59 noTheta // At end (because no pivot) 50 60 }; 51 61 /**@name Virtual method that the derived classes should provide. 
trunk/Clp/src/ClpNonLinearCost.cpp
r1723 r1761 1654 1654 } 1655 1655 if (CLP_METHOD2) { 1656 abort(); // may never work 1656 bound_[sequence]=0.0; 1657 cost2_[sequence]=costValue; 1658 setInitialStatus(status_[sequence]); 1657 1659 } 1658 1660 
trunk/Clp/src/ClpSimplex.cpp
r1738 r1761 8099 8099 upperIn_ = upper_[sequenceIn_]; 8100 8100 dualIn_ = dj_[sequenceIn_]; 8101 if (!nonLinearCost_) 8102 nonLinearCost_ = new ClpNonLinearCost(this); 8103 8101 8104 8102 8105 int returnCode = static_cast<ClpSimplexPrimal *> (this)>pivotResult(); … … 8108 8111 return 1; 8109 8112 } 8110 }8111 8112 /* Pivot out a variable and choose an incoing one. Assumes dual8113 feasible  will not go through a reduced cost.8114 Returns step length in theta8115 Returns ray in ray_ (or NULL if no pivot)8116 Return codes as before but 1 means no acceptable pivot8117 */8118 int8119 ClpSimplex::dualPivotResult()8120 {8121 return static_cast<ClpSimplexDual *> (this)>pivotResult();8122 8113 } 8123 8114 // Factorization frequency 
trunk/Clp/src/ClpSimplex.hpp
r1729 r1761 340 340 double * valueIncrease, int * sequenceIncrease, 341 341 double * valueDecrease, int * sequenceDecrease); 342 /** 343 Modifies coefficients etc and if necessary pivots in and out. 344 All at same status will be done (basis may go singular). 345 User can tell which others have been done (i.e. if status matches). 346 If called from outside will change status and return 0. 347 If called from event handler returns nonzero if user has to take action. 348 indices>=numberColumns are slacks (obviously no coefficients) 349 status array is (char) Status enum 350 */ 351 int modifyCoefficientsAndPivot(int number, 352 const int * which, 353 const CoinBigIndex * start, 354 const int * row, 355 const double * newCoefficient, 356 const unsigned char * newStatus=NULL, 357 const double * newLower=NULL, 358 const double * newUpper=NULL, 359 const double * newObjective=NULL); 342 360 /** Write the basis in MPS format to the specified file. 343 361 If writeValues true writes values of structurals … … 458 476 feasible  will not go through a reduced cost. 459 477 Returns step length in theta 460 Returns ray in ray_ (or NULL if no pivot)461 478 Return codes as before but 1 means no acceptable pivot 462 479 */ 463 int dualPivotResult(); 480 int dualPivotResultPart1(); 481 /** Do actual pivot 482 state is 0 if need tableau column, 1 if in rowArray_[1] 483 */ 484 int pivotResultPart2(int algorithm,int state); 464 485 465 486 /** Common bits of coding for dual and primal. Return 0 if okay, … … 839 860 return dualIn_; 840 861 } 862 /// Set reduced cost of last incoming to force error 863 inline void setDualIn(double value) { 864 dualIn_ = value; 865 } 841 866 /// Pivot Row for use by classes e.g. steepestedge 842 867 inline int pivotRow() const { … … 982 1007 valueOut_ = value; 983 1008 } 1009 /// Dual value of Out variable 1010 inline double dualOut() const { 1011 return dualOut_; 1012 } 1013 /// Set dual value of out variable 1014 inline void setDualOut(double value) { 1015 dualOut_ = value; 1016 } 984 1017 /// Set lower of out variable 985 1018 inline void setLowerOut(double value) { … … 1166 1199 inline int progressFlag() const { 1167 1200 return (progressFlag_ & 3); 1201 } 1202 /// For dealing with all issues of cycling etc 1203 inline ClpSimplexProgress * progress() 1204 { return &progress_;} 1205 /// Force refactorization early value 1206 inline int forceFactorization() const { 1207 return forceFactorization_ ; 1168 1208 } 1169 1209 /// Force refactorization early 
trunk/Clp/src/ClpSimplexDual.cpp
r1729 r1761 6443 6443 feasible  will not go through a reduced cost. 6444 6444 Returns step length in theta 6445 Returns ray in ray_ (or NULL if no pivot)6446 6445 Return codes as before but 1 means no acceptable pivot 6447 6446 */ 6448 6447 int 6449 ClpSimplexDual::pivotResult ()6448 ClpSimplexDual::pivotResultPart1() 6450 6449 { 6451 abort(); 6452 return 0; 6450 // Get good size for pivot 6451 // Allow first few iterations to take tiny 6452 double acceptablePivot = 1.0e1 * acceptablePivot_; 6453 if (numberIterations_ > 100) 6454 acceptablePivot = acceptablePivot_; 6455 if (factorization_>pivots() > 10) 6456 acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict 6457 else if (factorization_>pivots() > 5) 6458 acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict 6459 else if (factorization_>pivots()) 6460 acceptablePivot = acceptablePivot_; // relax 6461 // But factorizations complain if <1.0e8 6462 //acceptablePivot=CoinMax(acceptablePivot,1.0e8); 6463 double bestPossiblePivot = 1.0; 6464 // get sign for finding row of tableau 6465 // create as packed 6466 double direction = directionOut_; 6467 assert (!rowArray_[0]>getNumElements()); 6468 rowArray_[1]>clear(); //assert (!rowArray_[1]>getNumElements()); 6469 assert (!columnArray_[0]>getNumElements()); 6470 assert (!columnArray_[1]>getNumElements()); 6471 rowArray_[0]>createPacked(1, &pivotRow_, &direction); 6472 factorization_>updateColumnTranspose(rowArray_[1], rowArray_[0]); 6473 // Allow to do dualColumn0 6474 if (numberThreads_ < 1) 6475 spareIntArray_[0] = 1; 6476 spareDoubleArray_[0] = acceptablePivot; 6477 rowArray_[3]>clear(); 6478 sequenceIn_ = 1; 6479 // put row of tableau in rowArray[0] and columnArray[0] 6480 assert (!rowArray_[1]>getNumElements()); 6481 if (!scaledMatrix_) { 6482 if ((moreSpecialOptions_ & 8) != 0 && !rowScale_) 6483 spareIntArray_[0] = 1; 6484 matrix_>transposeTimes(this, 1.0, 6485 rowArray_[0], rowArray_[1], columnArray_[0]); 6486 } else { 6487 double * saveR = rowScale_; 6488 double * saveC = columnScale_; 6489 rowScale_ = NULL; 6490 columnScale_ = NULL; 6491 if ((moreSpecialOptions_ & 8) != 0) 6492 spareIntArray_[0] = 1; 6493 scaledMatrix_>transposeTimes(this, 1.0, 6494 rowArray_[0], rowArray_[1], columnArray_[0]); 6495 rowScale_ = saveR; 6496 columnScale_ = saveC; 6497 } 6498 // do ratio test for normal iteration 6499 dualOut_ *= 1.0e8; 6500 bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3], 6501 columnArray_[1], acceptablePivot, 6502 NULL/*dubiousWeights*/); 6503 dualOut_ *= 1.0e8; 6504 if (fabs(bestPossiblePivot)<1.0e6) 6505 return 1; 6506 else 6507 return 0; 6453 6508 } 6454 6509 /* 
trunk/Clp/src/ClpSimplexDual.hpp
r1665 r1761 282 282 /** Pivot in a variable and choose an outgoing one. Assumes dual 283 283 feasible  will not go through a reduced cost. Returns step length in theta 284 Returns ray in ray_ (or NULL if no pivot)285 284 Return codes as before but 1 means no acceptable pivot 286 285 */ 287 int pivotResult ();286 int pivotResultPart1(); 288 287 /** Get next free , 1 if none */ 289 288 int nextSuperBasic(); 
trunk/Clp/src/ClpSimplexOther.cpp
r1738 r1761 16 16 #include "ClpFactorization.hpp" 17 17 #include "ClpDualRowDantzig.hpp" 18 #include "ClpNonLinearCost.hpp" 18 19 #include "ClpDynamicMatrix.hpp" 19 20 #include "CoinPackedMatrix.hpp" … … 1889 1890 ClpDualRowPivot * savePivot = dualRowPivot_; 1890 1891 dualRowPivot_ = new ClpDualRowDantzig(); 1892 dualRowPivot_>setModel(this); 1891 1893 1892 1894 if (!returnCode) { … … 1983 1985 reinterpret_cast<ClpSimplexDual *> (this)>gutsOfDual(0, saveDuals, 1, data); 1984 1986 assert (!problemStatus_); 1987 for (int i=0;i<numberRows_+numberColumns_;i++) 1988 setFakeBound(i, noFake); 1985 1989 // Now do parametrics 1986 1990 handler_>message(CLP_PARAMETRICS_STATS, messages_) … … 2092 2096 << line << CoinMessageEol; 2093 2097 return problemStatus_; 2098 } 2099 /* Parametrics 2100 This is an initial slow version. 2101 The code uses current bounds + theta * change (if change array not NULL) 2102 It starts at startingTheta and returns ending theta in endingTheta. 2103 If it can not reach input endingTheta return code will be 1 for infeasible, 2104 2 for unbounded, if error on ranges 1, otherwise 0. 2105 Event handler may do more 2106 On exit endingTheta is maximum reached (can be used for next startingTheta) 2107 */ 2108 int 2109 ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, 2110 const double * changeLowerBound, const double * changeUpperBound, 2111 const double * changeLowerRhs, const double * changeUpperRhs) 2112 { 2113 int savePerturbation = perturbation_; 2114 perturbation_ = 102; // switch off 2115 algorithm_ = 1; 2116 2117 // save data 2118 ClpDataSave data = saveData(); 2119 // Dantzig (as will not be used) (out later) 2120 ClpDualRowPivot * savePivot = dualRowPivot_; 2121 dualRowPivot_ = new ClpDualRowDantzig(); 2122 dualRowPivot_>setModel(this); 2123 int numberTotal = numberRows_ + numberColumns_; 2124 /* 2125 Save information and modify 2126 */ 2127 double * saveLower = new double [3*numberTotal]; 2128 double * saveUpper = new double [3*numberTotal]; 2129 double * lowerCopy = saveLower+2*numberTotal; 2130 double * upperCopy = saveUpper+2*numberTotal; 2131 double * lowerChange = saveLower+numberTotal; 2132 double * upperChange = saveUpper+numberTotal; 2133 // Find theta when bounds will cross over and create arrays 2134 memset(lowerChange, 0, numberTotal * sizeof(double)); 2135 memset(upperChange, 0, numberTotal * sizeof(double)); 2136 if (changeLowerBound) 2137 memcpy(lowerChange,changeLowerBound,numberColumns_*sizeof(double)); 2138 if (changeUpperBound) 2139 memcpy(upperChange,changeUpperBound,numberColumns_*sizeof(double)); 2140 if (changeLowerRhs) 2141 memcpy(lowerChange+numberColumns_, 2142 changeLowerRhs,numberRows_*sizeof(double)); 2143 if (changeUpperRhs) 2144 memcpy(upperChange+numberColumns_, 2145 changeUpperRhs,numberRows_*sizeof(double)); 2146 memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double)); 2147 memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double)); 2148 memcpy(lowerCopy+numberColumns_, 2149 rowLower_,numberRows_*sizeof(double)); 2150 memcpy(upperCopy+numberColumns_, 2151 rowUpper_,numberRows_*sizeof(double)); 2152 assert (!rowScale_); 2153 double maxTheta = 1.0e50; 2154 for (int iRow = 0; iRow < numberRows_; iRow++) { 2155 double lower = rowLower_[iRow]; 2156 double upper = rowUpper_[iRow]; 2157 if (lower<1.0e30) 2158 lowerChange[numberColumns_+iRow]=0.0; 2159 double chgLower = lowerChange[numberColumns_+iRow]; 2160 if (upper>1.0e30) 2161 upperChange[numberColumns_+iRow]=0.0; 2162 double chgUpper = upperChange[numberColumns_+iRow]; 2163 if (lower > 1.0e30 && upper < 1.0e30) { 2164 if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) { 2165 maxTheta = (upper  lower) / (chgLower  chgUpper); 2166 } 2167 } 2168 lower+=startingTheta*chgLower; 2169 upper+=startingTheta*chgUpper; 2170 if (lower > upper) { 2171 maxTheta = 1.0; 2172 break; 2173 } 2174 rowLower_[iRow]=lower; 2175 rowUpper_[iRow]=upper; 2176 } 2177 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { 2178 double lower = columnLower_[iColumn]; 2179 double upper = columnUpper_[iColumn]; 2180 if (lower<1.0e30) 2181 lowerChange[iColumn]=0.0; 2182 double chgLower = lowerChange[iColumn]; 2183 if (upper>1.0e30) 2184 upperChange[iColumn]=0.0; 2185 double chgUpper = upperChange[iColumn]; 2186 if (lower > 1.0e30 && upper < 1.0e30) { 2187 if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) { 2188 maxTheta = (upper  lower) / (chgLower  chgUpper); 2189 } 2190 } 2191 lower+=startingTheta*chgLower; 2192 upper+=startingTheta*chgUpper; 2193 if (lower > upper) { 2194 maxTheta = 1.0; 2195 break; 2196 } 2197 columnLower_[iColumn]=lower; 2198 columnUpper_[iColumn]=upper; 2199 } 2200 if (maxTheta == 1.0e50) 2201 maxTheta = COIN_DBL_MAX; 2202 int returnCode=0; 2203 if (maxTheta < 0.0) { 2204 // bad ranges or initial 2205 returnCode = 1; 2206 } 2207 if (maxTheta < endingTheta) { 2208 char line[100]; 2209 sprintf(line,"Crossover considerations reduce ending theta from %g to %g\n", 2210 endingTheta,maxTheta); 2211 handler_>message(CLP_GENERAL,messages_) 2212 << line << CoinMessageEol; 2213 endingTheta = maxTheta; 2214 } 2215 if (endingTheta < startingTheta) { 2216 // bad initial 2217 returnCode = 2; 2218 } 2219 bool swapped=false; 2220 if (!returnCode) { 2221 returnCode = reinterpret_cast<ClpSimplexDual *> (this)>startupSolve(0, NULL, 0); 2222 if (!returnCode) { 2223 swapped=true; 2224 double * temp; 2225 memcpy(saveLower,lower_,numberTotal*sizeof(double)); 2226 temp=saveLower; 2227 saveLower=lower_; 2228 lower_=temp; 2229 memcpy(saveUpper,upper_,numberTotal*sizeof(double)); 2230 temp=saveUpper; 2231 saveUpper=upper_; 2232 upper_=temp; 2233 double saveEndingTheta = endingTheta; 2234 double * saveDuals = NULL; 2235 reinterpret_cast<ClpSimplexDual *> (this)>gutsOfDual(0, saveDuals, 1, data); 2236 assert (!problemStatus_); 2237 for (int i=0;i<numberRows_+numberColumns_;i++) 2238 setFakeBound(i, noFake); 2239 // Now do parametrics 2240 handler_>message(CLP_PARAMETRICS_STATS, messages_) 2241 << startingTheta << objectiveValue() << CoinMessageEol; 2242 while (!returnCode) { 2243 returnCode = parametricsLoop(startingTheta, endingTheta, 2244 data); 2245 if (!returnCode) { 2246 startingTheta = endingTheta; 2247 endingTheta = saveEndingTheta; 2248 handler_>message(CLP_PARAMETRICS_STATS, messages_) 2249 << startingTheta << objectiveValue() << CoinMessageEol; 2250 if (startingTheta >= endingTheta) 2251 break; 2252 } else if (returnCode == 1) { 2253 // trouble  do external solve 2254 abort(); //needToDoSomething = true; 2255 } else if (problemStatus_==1) { 2256 // can't move any further 2257 handler_>message(CLP_PARAMETRICS_STATS, messages_) 2258 << endingTheta << objectiveValue() << CoinMessageEol; 2259 problemStatus_=0; 2260 } 2261 } 2262 } 2263 reinterpret_cast<ClpSimplexDual *> (this)>finishSolve(0); 2264 if (swapped&&lower_) { 2265 double * temp=saveLower; 2266 saveLower=lower_; 2267 lower_=temp; 2268 temp=saveUpper; 2269 saveUpper=upper_; 2270 upper_=temp; 2271 } 2272 } 2273 memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double)); 2274 memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double)); 2275 memcpy(rowLower_,lowerCopy+numberColumns_, 2276 numberRows_*sizeof(double)); 2277 memcpy(rowUpper_,upperCopy+numberColumns_, 2278 numberRows_*sizeof(double)); 2279 delete [] saveLower; 2280 delete [] saveUpper; 2281 delete dualRowPivot_; 2282 dualRowPivot_ = savePivot; 2283 // Restore any saved stuff 2284 restoreData(data); 2285 perturbation_ = savePerturbation; 2286 char line[100]; 2287 sprintf(line,"Ending theta %g\n", endingTheta); 2288 handler_>message(CLP_GENERAL,messages_) 2289 << line << CoinMessageEol; 2290 return problemStatus_; 2094 2291 } 2095 2292 /* Version of parametrics which reads from file … … 2769 2966 } 2770 2967 } 2968 int 2969 ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, 2970 ClpDataSave & data) 2971 { 2972 int numberTotal = numberRows_+numberColumns_; 2973 const double * changeLower = lower_+numberTotal; 2974 const double * changeUpper = upper_+numberTotal; 2975 // stuff is already at starting 2976 // For this crude version just try and go to end 2977 double change = 0.0; 2978 int i; 2979 for ( i = 0; i < numberTotal; i++) { 2980 lower_[i] += change * changeLower[i]; 2981 upper_[i] += change * changeUpper[i]; 2982 switch(getStatus(i)) { 2983 2984 case basic: 2985 case isFree: 2986 case superBasic: 2987 break; 2988 case isFixed: 2989 case atUpperBound: 2990 solution_[i] = upper_[i]; 2991 break; 2992 case atLowerBound: 2993 solution_[i] = lower_[i]; 2994 break; 2995 } 2996 } 2997 problemStatus_ = 1; 2998 2999 // This says whether to restore things etc 3000 // startup will have factorized so can skip 3001 int factorType = 0; 3002 // Start check for cycles 3003 progress_.startCheck(); 3004 // Say change made on first iteration 3005 changeMade_ = 1; 3006 /* 3007 Status of problem: 3008 0  optimal 3009 1  infeasible 3010 2  unbounded 3011 1  iterating 3012 2  factorization wanted 3013 3  redo checking without factorization 3014 4  looks infeasible 3015 */ 3016 while (problemStatus_ < 0) { 3017 int iRow, iColumn; 3018 // clear 3019 for (iRow = 0; iRow < 4; iRow++) { 3020 rowArray_[iRow]>clear(); 3021 } 3022 3023 for (iColumn = 0; iColumn < 2; iColumn++) { 3024 columnArray_[iColumn]>clear(); 3025 } 3026 3027 // give matrix (and model costs and bounds a chance to be 3028 // refreshed (normally null) 3029 matrix_>refresh(this); 3030 // may factorize, checks if problem finished 3031 statusOfProblemInParametrics(factorType, data); 3032 // Say good factorization 3033 factorType = 1; 3034 if (data.sparseThreshold_) { 3035 // use default at present 3036 factorization_>sparseThreshold(0); 3037 factorization_>goSparse(); 3038 } 3039 3040 // exit if victory declared 3041 if (problemStatus_ >= 0 && startingTheta>=endingTheta1.0e7 ) 3042 break; 3043 3044 // test for maximum iterations 3045 if (hitMaximumIterations()) { 3046 problemStatus_ = 3; 3047 break; 3048 } 3049 // Check event 3050 { 3051 int status = eventHandler_>event(ClpEventHandler::endOfFactorization); 3052 if (status >= 0) { 3053 problemStatus_ = 5; 3054 secondaryStatus_ = ClpEventHandler::endOfFactorization; 3055 break; 3056 } 3057 } 3058 // Do iterations 3059 problemStatus_=1; 3060 whileIterating(startingTheta, endingTheta, 0.0, 3061 changeLower, changeUpper, 3062 NULL); 3063 startingTheta = endingTheta; 3064 } 3065 if (!problemStatus_) { 3066 theta_ = change + startingTheta; 3067 eventHandler_>event(ClpEventHandler::theta); 3068 return 0; 3069 } else if (problemStatus_ == 10) { 3070 return 1; 3071 } else { 3072 return problemStatus_; 3073 } 3074 } 2771 3075 /* Checks if finished. Updates status */ 2772 3076 void … … 2914 3218 // status 3 to go to top without an invert 2915 3219 int returnCode = 1; 2916 double saveSumDual = sumDualInfeasibilities_; // so we know to be careful2917 3220 double lastTheta = startingTheta; 2918 3221 double useTheta = startingTheta; … … 2929 3232 } 2930 3233 } 3234 #if 0 2931 3235 // See if objective 2932 3236 for (iSequence = 0; iSequence < numberTotal; iSequence++) { … … 2936 3240 } 2937 3241 } 2938 assert (type); 3242 #endif 3243 if (!type) { 3244 problemStatus_=0; 3245 return 0; 3246 } 2939 3247 while (problemStatus_ == 1) { 2940 3248 double increaseTheta = CoinMin(endingTheta  lastTheta, 1.0e50); … … 2946 3254 double change = useTheta  lastTheta; 2947 3255 for (int i = 0; i < numberTotal; i++) { 2948 lower_[i] += change * changeLower[i]; 2949 upper_[i] += change * changeUpper[i]; 2950 switch(getStatus(i)) { 2951 2952 case basic: 2953 case isFree: 2954 case superBasic: 2955 break; 2956 case isFixed: 2957 case atUpperBound: 2958 solution_[i] = upper_[i]; 2959 break; 2960 case atLowerBound: 2961 solution_[i] = lower_[i]; 2962 break; 3256 if (changeLower[i]changeUpper[i]) { 3257 lower_[i] += change * changeLower[i]; 3258 upper_[i] += change * changeUpper[i]; 3259 switch(getStatus(i)) { 3260 3261 case basic: 3262 case isFree: 3263 case superBasic: 3264 break; 3265 case isFixed: 3266 case atUpperBound: 3267 solution_[i] = upper_[i]; 3268 break; 3269 case atLowerBound: 3270 solution_[i] = lower_[i]; 3271 break; 3272 } 2963 3273 } 3274 #if 0 2964 3275 cost_[i] += change * changeObjective[i]; 3276 #endif 2965 3277 assert (solution_[i]>lower_[i]1.0e5&& 2966 3278 solution_[i]<upper_[i]+1.0e5); … … 2969 3281 if (pivotType) { 2970 3282 problemStatus_ = 2; 3283 if (!factorization_>pivots()&&pivotRow_<0) 3284 problemStatus_=2; 2971 3285 endingTheta = useTheta; 2972 3286 return 4; … … 2983 3297 // check accuracy of weights 2984 3298 dualRowPivot_>checkAccuracy(); 2985 // Get good size for pivot2986 // Allow first few iterations to take tiny2987 double acceptablePivot = 1.0e9;2988 if (numberIterations_ > 100)2989 acceptablePivot = 1.0e8;2990 if (factorization_>pivots() > 10 2991 (factorization_>pivots() && saveSumDual))2992 acceptablePivot = 1.0e5; // if we have iterated be more strict2993 else if (factorization_>pivots() > 5)2994 acceptablePivot = 1.0e6; // if we have iterated be slightly more strict2995 else if (factorization_>pivots())2996 acceptablePivot = 1.0e8; // relax2997 double bestPossiblePivot = 1.0;2998 // get sign for finding row of tableau2999 // normal iteration3000 // create as packed3001 double direction = directionOut_;3002 rowArray_[0]>createPacked(1, &pivotRow_, &direction);3003 factorization_>updateColumnTranspose(rowArray_[1], rowArray_[0]);3004 // put row of tableau in rowArray[0] and columnArray[0]3005 matrix_>transposeTimes(this, 1.0,3006 rowArray_[0], rowArray_[3], columnArray_[0]);3007 3299 // do ratio test for normal iteration 3008 bestPossiblePivot = reinterpret_cast<ClpSimplexDual *> ( this)>dualColumn(rowArray_[0], 3009 columnArray_[0], columnArray_[1], 3010 rowArray_[3], acceptablePivot, NULL); 3300 double bestPossiblePivot = bestPivot(); 3011 3301 if (sequenceIn_ >= 0) { 3012 3302 // normal iteration … … 3206 3496 objectiveChange = objectiveValue_; 3207 3497 // outgoing 3498 originalBound(sequenceOut_,useTheta,changeLower,changeUpper); 3499 lowerOut_=lower_[sequenceOut_]; 3500 upperOut_=upper_[sequenceOut_]; 3208 3501 // set dj to zero unless values pass 3209 3502 if (directionOut_ > 0) { … … 3254 3547 } 3255 3548 // and set bounds correctly 3256 reinterpret_cast<ClpSimplexDual *> ( this)>originalBound(sequenceIn_);3549 originalBound(sequenceIn_,useTheta,changeLower,changeUpper); 3257 3550 reinterpret_cast<ClpSimplexDual *> ( this)>changeBound(sequenceOut_); 3258 3551 if (whatNext == 1) { … … 3276 3569 } else { 3277 3570 // no incoming column is valid 3571 int savePivotRow = pivotRow_; 3278 3572 pivotRow_ = 1; 3279 3573 #ifdef CLP_DEBUG … … 3281 3575 printf("** no column pivot\n"); 3282 3576 #endif 3283 if (factorization_>pivots() < 5) { 3577 if (factorization_>pivots() < 5) { 3578 #if 0 3284 3579 // If not in branch and bound etc save ray 3285 3580 if ((specialOptions_&(1024  4096)) == 0) { … … 3290 3585 ClpDisjointCopyN(rowArray_[0]>denseVector(), numberRows_, ray_); 3291 3586 } 3587 #endif 3292 3588 // If we have just factorized and infeasibility reasonable say infeas 3293 3589 if (((specialOptions_ & 4096) != 0  bestPossiblePivot < 1.0e11) && dualBound_ > 1.0e8) { … … 3303 3599 rowArray_[0]>clear(); 3304 3600 columnArray_[0]>clear(); 3601 pivotRow_ = savePivotRow; 3602 int action = eventHandler_>event(ClpEventHandler::noTheta); 3603 if (action==2) 3604 problemStatus_=1; // carry on 3305 3605 returnCode = 1; 3306 3606 break; … … 3453 3753 delete [] primalChange; 3454 3754 delete [] dualChange; 3455 endingTheta = lastTheta ;3755 endingTheta = lastTheta+theta_; 3456 3756 return returnCode; 3757 } 3758 // Finds best possible pivot 3759 double 3760 ClpSimplexOther::bestPivot(bool justColumns) 3761 { 3762 // Get good size for pivot 3763 // Allow first few iterations to take tiny 3764 double acceptablePivot = 1.0e9; 3765 if (numberIterations_ > 100) 3766 acceptablePivot = 1.0e8; 3767 if (factorization_>pivots() > 10  3768 (factorization_>pivots() && sumDualInfeasibilities_)) 3769 acceptablePivot = 1.0e5; // if we have iterated be more strict 3770 else if (factorization_>pivots() > 5) 3771 acceptablePivot = 1.0e6; // if we have iterated be slightly more strict 3772 else if (factorization_>pivots()) 3773 acceptablePivot = 1.0e8; // relax 3774 double bestPossiblePivot = 1.0; 3775 // get sign for finding row of tableau 3776 // normal iteration 3777 // create as packed 3778 double direction = directionOut_; 3779 rowArray_[0]>createPacked(1, &pivotRow_, &direction); 3780 factorization_>updateColumnTranspose(rowArray_[1], rowArray_[0]); 3781 // put row of tableau in rowArray[0] and columnArray[0] 3782 matrix_>transposeTimes(this, 1.0, 3783 rowArray_[0], rowArray_[3], columnArray_[0]); 3784 sequenceIn_=1; 3785 if (justColumns) 3786 rowArray_[0]>clear(); 3787 // do ratio test for normal iteration 3788 bestPossiblePivot = 3789 reinterpret_cast<ClpSimplexDual *> 3790 ( this)>dualColumn(rowArray_[0], 3791 columnArray_[0], columnArray_[1], 3792 rowArray_[3], acceptablePivot, NULL); 3793 return bestPossiblePivot; 3457 3794 } 3458 3795 // Computes next theta and says if objective or bounds (0= bounds, 1 objective, 1 none) … … 3569 3906 return 0; 3570 3907 } else { 3908 theta_=0.0; 3571 3909 return 1; 3910 } 3911 } 3912 // Restores bound to original bound 3913 void 3914 ClpSimplexOther::originalBound(int iSequence, double theta, 3915 const double * changeLower, 3916 const double * changeUpper) 3917 { 3918 if (getFakeBound(iSequence) != noFake) { 3919 numberFake_; 3920 setFakeBound(iSequence, noFake); 3921 if (iSequence >= numberColumns_) { 3922 // rows 3923 int iRow = iSequence  numberColumns_; 3924 rowLowerWork_[iRow] = rowLower_[iRow]+theta*changeLower[iSequence]; 3925 rowUpperWork_[iRow] = rowUpper_[iRow]+theta*changeUpper[iSequence]; 3926 if (rowScale_) { 3927 if (rowLowerWork_[iRow] > 1.0e50) 3928 rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_; 3929 if (rowUpperWork_[iRow] < 1.0e50) 3930 rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_; 3931 } else if (rhsScale_ != 1.0) { 3932 if (rowLowerWork_[iRow] > 1.0e50) 3933 rowLowerWork_[iRow] *= rhsScale_; 3934 if (rowUpperWork_[iRow] < 1.0e50) 3935 rowUpperWork_[iRow] *= rhsScale_; 3936 } 3937 } else { 3938 // columns 3939 columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*changeLower[iSequence]; 3940 columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*changeUpper[iSequence]; 3941 if (rowScale_) { 3942 double multiplier = 1.0 * inverseColumnScale_[iSequence]; 3943 if (columnLowerWork_[iSequence] > 1.0e50) 3944 columnLowerWork_[iSequence] *= multiplier * rhsScale_; 3945 if (columnUpperWork_[iSequence] < 1.0e50) 3946 columnUpperWork_[iSequence] *= multiplier * rhsScale_; 3947 } else if (rhsScale_ != 1.0) { 3948 if (columnLowerWork_[iSequence] > 1.0e50) 3949 columnLowerWork_[iSequence] *= rhsScale_; 3950 if (columnUpperWork_[iSequence] < 1.0e50) 3951 columnUpperWork_[iSequence] *= rhsScale_; 3952 } 3953 } 3572 3954 } 3573 3955 } … … 5006 5388 //printf("objective value is %g\n", objValue); 5007 5389 } 5390 /* 5391 Modifies coefficients etc and if necessary pivots in and out. 5392 All at same status will be done (basis may go singular). 5393 User can tell which others have been done (i.e. if status matches). 5394 If called from outside will change status and return 0 (1 if not ClpPacked) 5395 If called from event handler returns nonzero if user has to take action. 5396 1  not ClpPackedMatrix 5397 0  no pivots 5398 1  pivoted 5399 2  pivoted and optimal 5400 3  refactorize 5401 4  worse than that 5402 indices>=numberColumns are slacks (obviously no coefficients) 5403 status array is (char) Status enum 5404 */ 5405 int 5406 ClpSimplex::modifyCoefficientsAndPivot(int number, 5407 const int * which, 5408 const CoinBigIndex * start, 5409 const int * row, 5410 const double * newCoefficient, 5411 const unsigned char * newStatus, 5412 const double * newLower, 5413 const double * newUpper, 5414 const double * newObjective) 5415 { 5416 ClpPackedMatrix* clpMatrix = 5417 dynamic_cast< ClpPackedMatrix*>(matrix_); 5418 bool canPivot = lower_!=NULL && factorization_!=NULL; 5419 int returnCode=0; 5420 if (!clpMatrix) { 5421 canPivot=false; 5422 returnCode=1; 5423 // very slow 5424 for (int i=0;i<number;i++) { 5425 int iSequence=which[i]; 5426 if (iSequence<numberColumns_) { 5427 for (CoinBigIndex j=start[i];j<start[i+1];j++) { 5428 matrix_>modifyCoefficient(row[j],iSequence,newCoefficient[j]); 5429 } 5430 } else { 5431 assert (start[i]==start[i+1]); 5432 } 5433 } 5434 } else { 5435 CoinPackedMatrix * matrix = clpMatrix>getPackedMatrix(); 5436 matrix>modifyCoefficients(number,which,start, 5437 row,newCoefficient); 5438 if (canPivot) { 5439 // ? faster to modify row copy?? 5440 if (rowCopy_&&start[number]) { 5441 delete rowCopy_; 5442 rowCopy_ = clpMatrix>reverseOrderedCopy(); 5443 } 5444 assert (!newStatus); // do later 5445 int numberPivots = factorization_>pivots(); 5446 int needed=0; 5447 for (int i=0;i<number;i++) { 5448 int iSequence=which[i]; 5449 if (start[i+1]>start[i]&&getStatus(iSequence)==basic) 5450 needed++; 5451 } 5452 if (needed&&numberPivots+needed<20&&needed<2) { 5453 // update factorization 5454 int saveIn = sequenceIn_; 5455 int savePivot = pivotRow_; 5456 int nArray=0; 5457 CoinIndexedVector * vec[2]; 5458 for (int i=0;i<4;i++) { 5459 if (!rowArray_[i]>getNumElements()) { 5460 vec[nArray++]=rowArray_[i]; 5461 if (nArray==2) 5462 break; 5463 } 5464 } 5465 assert (nArray==2); // could use temp array 5466 for (int i=0;i<number;i++) { 5467 int sequenceIn_=which[i]; 5468 if (start[i+1]>start[i]&&getStatus(sequenceIn_)==basic) { 5469 // find pivot row 5470 for (pivotRow_=0;pivotRow_<numberRows_;pivotRow_++) { 5471 if (pivotVariable_[pivotRow_]==sequenceIn_) 5472 break; 5473 } 5474 assert(pivotRow_<numberRows_); 5475 // unpack column 5476 assert(!vec[0]>getNumElements()); 5477 unpackPacked(vec[0]); 5478 // update 5479 assert(!vec[1]>getNumElements()); 5480 factorization_>updateColumnFT(vec[1], vec[0]); 5481 // Find alpha 5482 const int * indices = vec[0]>getIndices(); 5483 const double * array = vec[0]>denseVector(); 5484 int n=vec[0]>getNumElements(); 5485 alpha_=0.0; 5486 for (int i=0;i<n;i++) { 5487 if (indices[i]==pivotRow_) { 5488 alpha_ = array[i]; 5489 break; 5490 } 5491 } 5492 int updateStatus=2; 5493 if (fabs(alpha_)>1.0e7) 5494 updateStatus = factorization_>replaceColumn(this, 5495 vec[1], 5496 vec[0], 5497 pivotRow_, 5498 alpha_); 5499 assert(!vec[1]>getNumElements()); 5500 vec[0]>clear(); 5501 if (updateStatus) { 5502 returnCode=3; 5503 break; 5504 } 5505 } 5506 } 5507 sequenceIn_=saveIn; 5508 pivotRow_ = savePivot; 5509 if (!returnCode) 5510 returnCode=100; // say can do more 5511 } else if (needed) { 5512 returnCode=3; // refactorize 5513 } 5514 } 5515 } 5516 if (newStatus) { 5517 for (int i=0;i<number;i++) { 5518 int iSequence=which[i]; 5519 status_[iSequence]=newStatus[i]; 5520 } 5521 } 5522 if (newLower) { 5523 for (int i=0;i<number;i++) { 5524 int iSequence=which[i]; 5525 if (iSequence<numberColumns_) { 5526 if (columnLower_[iSequence]!=newLower[i]) { 5527 columnLower_[iSequence]=newLower[i]; 5528 } 5529 } else { 5530 iSequence = numberColumns_; 5531 if (rowLower_[iSequence]!=newLower[i]) { 5532 rowLower_[iSequence]=newLower[i]; 5533 } 5534 } 5535 } 5536 } 5537 if (newLower) { 5538 for (int i=0;i<number;i++) { 5539 int iSequence=which[i]; 5540 if (iSequence<numberColumns_) { 5541 if (columnLower_[iSequence]!=newLower[i]) { 5542 columnLower_[iSequence]=newLower[i]; 5543 } 5544 } else { 5545 iSequence = numberColumns_; 5546 if (rowLower_[iSequence]!=newLower[i]) { 5547 rowLower_[iSequence]=newLower[i]; 5548 } 5549 } 5550 } 5551 } 5552 if (newUpper) { 5553 for (int i=0;i<number;i++) { 5554 int iSequence=which[i]; 5555 if (iSequence<numberColumns_) { 5556 if (columnUpper_[iSequence]!=newUpper[i]) { 5557 columnUpper_[iSequence]=newUpper[i]; 5558 } 5559 } else { 5560 iSequence = numberColumns_; 5561 if (rowUpper_[iSequence]!=newUpper[i]) { 5562 rowUpper_[iSequence]=newUpper[i]; 5563 } 5564 } 5565 } 5566 } 5567 if (newObjective) { 5568 double * obj = objective(); 5569 for (int i=0;i<number;i++) { 5570 int iSequence=which[i]; 5571 if (iSequence<numberColumns_) { 5572 if (obj[iSequence]!=newObjective[i]) { 5573 obj[iSequence]=newObjective[i]; 5574 } 5575 } else { 5576 assert (!newObjective[i]); 5577 } 5578 } 5579 } 5580 if (canPivot) { 5581 // update lower, upper, objective and nonLinearCost 5582 assert (!rowScale_); // for now 5583 memcpy(lower_,columnLower_,numberColumns_*sizeof(double)); 5584 memcpy(lower_+numberColumns_,rowLower_,numberRows_*sizeof(double)); 5585 memcpy(upper_,columnUpper_,numberColumns_*sizeof(double)); 5586 memcpy(upper_+numberColumns_,rowUpper_,numberRows_*sizeof(double)); 5587 memcpy(cost_,objective(),numberColumns_*sizeof(double)); 5588 memset(cost_+numberColumns_,0,numberRows_*sizeof(double)); 5589 // ? parameter to say no gutsOfSolution needed 5590 // make sure slacks have correct value 5591 // see if still optimal 5592 if (returnCode==100) { 5593 // is this needed 5594 if (nonLinearCost_) { 5595 // speed up later 5596 //nonLinearCost_>checkInfeasibilities(oldTolerance); 5597 delete nonLinearCost_; 5598 nonLinearCost_ = new ClpNonLinearCost(this); 5599 } 5600 gutsOfSolution(NULL,NULL,false); 5601 assert (!newStatus); 5602 printf("%d primal %d dual\n",numberPrimalInfeasibilities_, 5603 numberDualInfeasibilities_); 5604 returnCode=3; 5605 } else { 5606 // is this needed 5607 if (nonLinearCost_) { 5608 // speed up later 5609 #if 1 5610 for (int i=0;i<number;i++) { 5611 int iSequence=which[i]; 5612 nonLinearCost_>setOne(iSequence,solution_[iSequence], 5613 lower_[iSequence],upper_[iSequence], 5614 cost_[iSequence]); 5615 } 5616 #else 5617 //nonLinearCost_>checkInfeasibilities(oldTolerance); 5618 delete nonLinearCost_; 5619 nonLinearCost_ = new ClpNonLinearCost(this); 5620 //nonLinearCost_>checkInfeasibilities(0.0); 5621 #endif 5622 //gutsOfSolution(NULL,NULL,false); 5623 assert (!newStatus); 5624 } 5625 } 5626 } 5627 return returnCode; 5628 } 5629 5630 /* Pivot out a variable and choose an incoing one. Assumes dual 5631 feasible  will not go through a reduced cost. 5632 Returns step length in theta 5633 Return codes as before but 1 means no acceptable pivot 5634 */ 5635 int 5636 ClpSimplex::dualPivotResultPart1() 5637 { 5638 return static_cast<ClpSimplexDual *> (this)>pivotResultPart1(); 5639 } 5640 /* Do actual pivot 5641 state is 1,3 if got tableau column in rowArray_[1] 5642 2,3 if got tableau row in rowArray_[0] and columnArray_[0] 5643 */ 5644 int 5645 ClpSimplex::pivotResultPart2(int algorithm,int state) 5646 { 5647 if (!(state&1)) { 5648 // update the incoming column 5649 unpackPacked(rowArray_[1]); 5650 factorization_>updateColumnFT(rowArray_[2], rowArray_[1]); 5651 } 5652 #define CHECK_TABLEAU 0 5653 if (!(state&2)CHECK_TABLEAU) { 5654 // get tableau row 5655 // create as packed 5656 double direction = directionOut_; 5657 assert (!rowArray_[2]>getNumElements()); 5658 assert (!columnArray_[1]>getNumElements()); 5659 #if CHECK_TABLEAU 5660 printf("rowArray0 old\n"); 5661 rowArray_[0]>print(); 5662 rowArray_[0]>clear(); 5663 printf("columnArray0 old\n"); 5664 columnArray_[0]>print(); 5665 columnArray_[0]>clear(); 5666 #else 5667 assert (!columnArray_[0]>getNumElements()); 5668 assert (!rowArray_[0]>getNumElements()); 5669 #endif 5670 rowArray_[0]>createPacked(1, &pivotRow_, &direction); 5671 factorization_>updateColumnTranspose(rowArray_[2], rowArray_[0]); 5672 rowArray_[3]>clear(); 5673 // put row of tableau in rowArray[0] and columnArray[0] 5674 assert (!rowArray_[2]>getNumElements()); 5675 matrix_>transposeTimes(this, 1.0, 5676 rowArray_[0], rowArray_[2], columnArray_[0]); 5677 #if CHECK_TABLEAU 5678 printf("rowArray0 new\n"); 5679 rowArray_[0]>print(); 5680 printf("columnArray0 new\n"); 5681 columnArray_[0]>print(); 5682 #endif 5683 } 5684 assert (pivotRow_>=0); 5685 assert (sequenceIn_>=0); 5686 assert (sequenceOut_>=0); 5687 int returnCode=1; 5688 if (algorithm>0) { 5689 // replace in basis 5690 int updateStatus = factorization_>replaceColumn(this, 5691 rowArray_[2], 5692 rowArray_[1], 5693 pivotRow_, 5694 alpha_); 5695 if (!updateStatus) { 5696 dualIn_ = cost_[sequenceIn_]; 5697 double * work = rowArray_[1]>denseVector(); 5698 int number = rowArray_[1]>getNumElements(); 5699 int * which = rowArray_[1]>getIndices(); 5700 for (int i = 0; i < number; i++) { 5701 int iRow = which[i]; 5702 double alpha = work[i]; 5703 int iPivot = pivotVariable_[iRow]; 5704 dualIn_ = alpha * cost_[iPivot]; 5705 } 5706 returnCode=0; 5707 double multiplier = dualIn_ / alpha_; 5708 // update column djs 5709 int i; 5710 int * index = columnArray_[0]>getIndices(); 5711 number = columnArray_[0]>getNumElements(); 5712 double * element = columnArray_[0]>denseVector(); 5713 assert (columnArray_[0]>packedMode()); 5714 for (i = 0; i < number; i++) { 5715 int iSequence = index[i]; 5716 dj_[iSequence] += multiplier*element[i]; 5717 reducedCost_[iSequence] = dj_[iSequence]; 5718 element[i] = 0.0; 5719 } 5720 columnArray_[0]>setNumElements(0); 5721 // and row djs 5722 index = rowArray_[0]>getIndices(); 5723 number = rowArray_[0]>getNumElements(); 5724 element = rowArray_[0]>denseVector(); 5725 assert (rowArray_[0]>packedMode()); 5726 for (i = 0; i < number; i++) { 5727 int iSequence = index[i]; 5728 dj_[iSequence+numberColumns_] += multiplier*element[i]; 5729 dual_[iSequence] = dj_[iSequence+numberColumns_]; 5730 element[i] = 0.0; 5731 } 5732 rowArray_[0]>setNumElements(0); 5733 double oldCost = cost_[sequenceOut_]; 5734 // update primal solution 5735 5736 double objectiveChange = 0.0; 5737 // after this rowArray_[1] is not empty  used to update djs 5738 static_cast<ClpSimplexPrimal *>(this)>updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 0); 5739 double oldValue = valueIn_; 5740 if (directionIn_ == 1) { 5741 // as if from upper bound 5742 if (sequenceIn_ != sequenceOut_) { 5743 // variable becoming basic 5744 valueIn_ = fabs(theta_); 5745 } else { 5746 valueIn_ = lowerIn_; 5747 } 5748 } else { 5749 // as if from lower bound 5750 if (sequenceIn_ != sequenceOut_) { 5751 // variable becoming basic 5752 valueIn_ += fabs(theta_); 5753 } else { 5754 valueIn_ = upperIn_; 5755 } 5756 } 5757 objectiveChange += dualIn_ * (valueIn_  oldValue); 5758 // outgoing 5759 if (sequenceIn_ != sequenceOut_) { 5760 if (directionOut_ > 0) { 5761 valueOut_ = lowerOut_; 5762 } else { 5763 valueOut_ = upperOut_; 5764 } 5765 if(valueOut_ < lower_[sequenceOut_]  primalTolerance_) 5766 valueOut_ = lower_[sequenceOut_]  0.9 * primalTolerance_; 5767 else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_) 5768 valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_; 5769 // may not be exactly at bound and bounds may have changed 5770 // Make sure outgoing looks feasible 5771 directionOut_ = nonLinearCost_>setOneOutgoing(sequenceOut_, valueOut_); 5772 // May have got inaccurate 5773 //if (oldCost!=cost_[sequenceOut_]) 5774 //printf("costchange on %d from %g to %g\n",sequenceOut_, 5775 // oldCost,cost_[sequenceOut_]); 5776 dj_[sequenceOut_] = cost_[sequenceOut_]  oldCost; // normally updated next iteration 5777 solution_[sequenceOut_] = valueOut_; 5778 } 5779 // change cost and bounds on incoming if primal 5780 nonLinearCost_>setOne(sequenceIn_, valueIn_); 5781 progress_.startCheck(); // make sure won't worry about cycling 5782 int whatNext = housekeeping(objectiveChange); 5783 if (whatNext == 1) { 5784 returnCode = 2; // refactorize 5785 } else if (whatNext == 2) { 5786 // maximum iterations or equivalent 5787 returnCode = 3; 5788 } else if(numberIterations_ == lastGoodIteration_ 5789 + 2 * factorization_>maximumPivots()) { 5790 // done a lot of flips  be safe 5791 returnCode = 2; // refactorize 5792 } 5793 } else { 5794 // ? 5795 abort(); 5796 } 5797 } else { 5798 // dual 5799 // recompute dualOut_ 5800 if (directionOut_ < 0) { 5801 dualOut_ = valueOut_  upperOut_; 5802 } else { 5803 dualOut_ = lowerOut_  valueOut_; 5804 } 5805 // update the incoming column 5806 double btranAlpha = alpha_ * directionOut_; // for check 5807 rowArray_[1]>clear(); 5808 unpackPacked(rowArray_[1]); 5809 // moved into updateWeights  factorization_>updateColumnFT(rowArray_[2],rowArray_[1]); 5810 // and update dual weights (can do in parallel  with extra array) 5811 alpha_ = dualRowPivot_>updateWeights(rowArray_[0], 5812 rowArray_[2], 5813 rowArray_[3], 5814 rowArray_[1]); 5815 // see if update stable 5816 #ifdef CLP_DEBUG 5817 if ((handler_>logLevel() & 32)) 5818 printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_); 5819 #endif 5820 double checkValue = 1.0e7; 5821 // if can't trust much and long way from optimal then relax 5822 if (largestPrimalError_ > 10.0) 5823 checkValue = CoinMin(1.0e4, 1.0e8 * largestPrimalError_); 5824 if (fabs(btranAlpha) < 1.0e12  fabs(alpha_) < 1.0e12  5825 fabs(btranAlpha  alpha_) > checkValue*(1.0 + fabs(alpha_))) { 5826 handler_>message(CLP_DUAL_CHECK, messages_) 5827 << btranAlpha 5828 << alpha_ 5829 << CoinMessageEol; 5830 if (factorization_>pivots()) { 5831 dualRowPivot_>unrollWeights(); 5832 problemStatus_ = 2; // factorize now 5833 rowArray_[0]>clear(); 5834 rowArray_[1]>clear(); 5835 columnArray_[0]>clear(); 5836 returnCode = 2; 5837 abort(); 5838 return returnCode; 5839 } else { 5840 // take on more relaxed criterion 5841 double test; 5842 if (fabs(btranAlpha) < 1.0e8  fabs(alpha_) < 1.0e8) 5843 test = 1.0e1 * fabs(alpha_); 5844 else 5845 test = 1.0e4 * (1.0 + fabs(alpha_)); 5846 if (fabs(btranAlpha) < 1.0e12  fabs(alpha_) < 1.0e12  5847 fabs(btranAlpha  alpha_) > test) { 5848 abort(); 5849 } 5850 } 5851 } 5852 // update duals BEFORE replaceColumn so can do updateColumn 5853 double objectiveChange = 0.0; 5854 // do duals first as variables may flip bounds 5855 // rowArray_[0] and columnArray_[0] may have flips 5856 // so use rowArray_[3] for work array from here on 5857 int nswapped = 0; 5858 //rowArray_[0]>cleanAndPackSafe(1.0e60); 5859 //columnArray_[0]>cleanAndPackSafe(1.0e60); 5860 // make sure incoming doesn't count 5861 Status saveStatus = getStatus(sequenceIn_); 5862 setStatus(sequenceIn_, basic); 5863 nswapped = 5864 static_cast<ClpSimplexDual *>(this)>updateDualsInDual(rowArray_[0], columnArray_[0], 5865 rowArray_[2], theta_, 5866 objectiveChange, false); 5867 assert (!nswapped); 5868 setStatus(sequenceIn_, saveStatus); 5869 double oldDualOut = dualOut_; 5870 // which will change basic solution 5871 if (nswapped) { 5872 if (rowArray_[2]>getNumElements()) { 5873 factorization_>updateColumn(rowArray_[3], rowArray_[2]); 5874 dualRowPivot_>updatePrimalSolution(rowArray_[2], 5875 1.0, objectiveChange); 5876 } 5877 // recompute dualOut_ 5878 valueOut_ = solution_[sequenceOut_]; 5879 if (directionOut_ < 0) { 5880 dualOut_ = valueOut_  upperOut_; 5881 } else { 5882 dualOut_ = lowerOut_  valueOut_; 5883 } 5884 } 5885 // amount primal will move 5886 double movement = dualOut_ * directionOut_ / alpha_; 5887 double movementOld = oldDualOut * directionOut_ / alpha_; 5888 // so objective should increase by fabs(dj)*movement 5889 // but we already have objective change  so check will be good 5890 if (objectiveChange + fabs(movementOld * dualIn_) < CoinMax(1.0e5, 1.0e12 * fabs(objectiveValue_))) { 5891 if (handler_>logLevel() & 32) 5892 printf("movement %g, swap change %g, rest %g * %g\n", 5893 objectiveChange + fabs(movement * dualIn_), 5894 objectiveChange, movement, dualIn_); 5895 } 5896 // if stable replace in basis 5897 int updateStatus = factorization_>replaceColumn(this, 5898 rowArray_[2], 5899 rowArray_[1], 5900 pivotRow_, 5901 alpha_); 5902 // If looks like bad pivot  refactorize 5903 if (fabs(dualOut_) > 1.0e50) 5904 updateStatus = 2; 5905 // if no pivots, bad update but reasonable alpha  take and invert 5906 if (updateStatus == 2 && 5907 !factorization_>pivots() && fabs(alpha_) > 1.0e5) 5908 updateStatus = 4; 5909 if (updateStatus == 1  updateStatus == 4) { 5910 // slight error 5911 if (factorization_>pivots() > 5  updateStatus == 4) { 5912 problemStatus_ = 2; // factorize now 5913 returnCode = 3; 5914 } 5915 } else if (updateStatus == 2) { 5916 // major error 5917 dualRowPivot_>unrollWeights(); 5918 // later we may need to unwind more e.g. fake bounds 5919 if (factorization_>pivots() && 5920 ((moreSpecialOptions_ & 16) == 0  factorization_>pivots() > 4)) { 5921 problemStatus_ = 2; // factorize now 5922 returnCode = 2; 5923 moreSpecialOptions_ = 16; 5924 return returnCode; 5925 } else { 5926 // need to reject something 5927 abort(); 5928 } 5929 } else if (updateStatus == 3) { 5930 // out of memory 5931 // increase space if not many iterations 5932 if (factorization_>pivots() < 5933 0.5 * factorization_>maximumPivots() && 5934 factorization_>pivots() < 200) 5935 factorization_>areaFactor( 5936 factorization_>areaFactor() * 1.1); 5937 problemStatus_ = 2; // factorize now 5938 } else if (updateStatus == 5) { 5939 problemStatus_ = 2; // factorize now 5940 } 5941 // update primal solution 5942 if (theta_ < 0.0) { 5943 if (handler_>logLevel() & 32) 5944 printf("negative theta %g\n", theta_); 5945 theta_ = 0.0; 5946 } 5947 // do actual flips (should not be any?) 5948 static_cast<ClpSimplexDual *>(this)>flipBounds(rowArray_[0], columnArray_[0]); 5949 //rowArray_[1]>expand(); 5950 dualRowPivot_>updatePrimalSolution(rowArray_[1], 5951 movement, 5952 objectiveChange); 5953 // modify dualout 5954 dualOut_ /= alpha_; 5955 dualOut_ *= directionOut_; 5956 //setStatus(sequenceIn_,basic); 5957 dj_[sequenceIn_] = 0.0; 5958 double oldValue = valueIn_; 5959 if (directionIn_ == 1) { 5960 // as if from upper bound 5961 valueIn_ = upperIn_ + dualOut_; 5962 } else { 5963 // as if from lower bound 5964 valueIn_ = lowerIn_ + dualOut_; 5965 } 5966 objectiveChange += cost_[sequenceIn_] * (valueIn_  oldValue); 5967 // outgoing 5968 // set dj to zero unless values pass 5969 if (directionOut_ > 0) { 5970 valueOut_ = lowerOut_; 5971 dj_[sequenceOut_] = theta_; 5972 } else { 5973 valueOut_ = upperOut_; 5974 dj_[sequenceOut_] = theta_; 5975 } 5976 solution_[sequenceOut_] = valueOut_; 5977 int whatNext = housekeeping(objectiveChange); 5978 // and set bounds correctly 5979 static_cast<ClpSimplexDual *>(this)>originalBound(sequenceIn_); 5980 static_cast<ClpSimplexDual *>(this)>changeBound(sequenceOut_); 5981 if (whatNext == 1) { 5982 problemStatus_ = 2; // refactorize 5983 } else if (whatNext == 2) { 5984 // maximum iterations or equivalent 5985 problemStatus_ = 3; 5986 returnCode = 3; 5987 abort(); 5988 } 5989 } 5990 // Check event 5991 { 5992 int status = eventHandler_>event(ClpEventHandler::endOfIteration); 5993 if (status >= 0) { 5994 problemStatus_ = 5; 5995 secondaryStatus_ = ClpEventHandler::endOfIteration; 5996 returnCode = 3; 5997 } 5998 } 5999 // need to be able to refactoriza 6000 //printf("return code %d problem status %d\n", 6001 // returnCode,problemStatus_); 6002 return returnCode; 6003 } 
trunk/Clp/src/ClpSimplexOther.hpp
r1678 r1761 89 89 Returns 2 if unable to open file */ 90 90 int parametrics(const char * dataFile); 91 /** Parametrics 92 This is an initial slow version. 93 The code uses current bounds + theta * change (if change array not NULL) 94 It starts at startingTheta and returns ending theta in endingTheta. 95 If it can not reach input endingTheta return code will be 1 for infeasible, 96 2 for unbounded, if error on ranges 1, otherwise 0. 97 Event handler may do more 98 On exit endingTheta is maximum reached (can be used for next startingTheta) 99 */ 100 int parametrics(double startingTheta, double & endingTheta, 101 const double * changeLowerBound, const double * changeUpperBound, 102 const double * changeLowerRhs, const double * changeUpperRhs); 103 /// Finds best possible pivot 104 double bestPivot(bool justColumns=false); 91 105 92 106 private: … … 103 117 const double * changeObjective, ClpDataSave & data, 104 118 bool canTryQuick); 119 int parametricsLoop(double startingTheta, double & endingTheta, 120 ClpDataSave & data); 105 121 /** Refactorizes if necessary 106 122 Checks if finished. Updates status. … … 131 147 const double * changeLower, const double * changeUpper, 132 148 const double * changeObjective); 149 /// Restores bound to original bound 150 void originalBound(int iSequence, double theta, const double * changeLower, 151 const double * changeUpper); 133 152 /** 134 153 Row array has row part of pivot row 
trunk/Clp/src/ClpSimplexPrimal.cpp
r1738 r1761 460 460 461 461 // exit if victory declared 462 if (problemStatus_ >= 0) 463 break; 462 if (problemStatus_ >= 0) { 463 #ifdef CLP_USER_DRIVEN 464 int status = 465 eventHandler_>event(ClpEventHandler::endInPrimal); 466 if (status>=0&&status<10) { 467 // carry on 468 problemStatus_=1; 469 if (status==0) 470 continue; // refactorize 471 } else if (status>=10) { 472 problemStatus_=status10; 473 break; 474 } else { 475 break; 476 } 477 #else 478 break; 479 #endif 480 } 464 481 465 482 // test for maximum iterations … … 720 737 // Force to refactorize early next time 721 738 int numberPivots = factorization_>pivots(); 722 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); 723 returnCode = 0; 724 break; 739 returnCode = 0; 740 #ifdef CLP_USER_DRIVEN 741 // If large number of pivots trap later? 742 if (problemStatus_==1 && numberPivots<1000) { 743 int status = eventHandler_>event(ClpEventHandler::noCandidateInPrimal); 744 if (status>=0&&status<10) { 745 // carry on 746 problemStatus_=1; 747 if (status==0) 748 break; 749 } else if (status>=10) { 750 problemStatus_=status10; 751 break; 752 } else { 753 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); 754 break; 755 } 756 } 757 #else 758 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1); 759 break; 760 #endif 725 761 } 726 762 } … … 856 892 static_cast<double>(numberDualInfeasibilities_ + 10); 857 893 #endif 894 #ifdef CLP_USER_DRIVEN 895 int status = eventHandler_>event(ClpEventHandler::goodFactorization); 896 if (status >= 0) { 897 numberThrownOut=status; 898 } else { 899 numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0)); 900 } 901 #else 858 902 numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0)); 903 #endif 859 904 double sumInfeasibility = nonLinearCost_>sumInfeasibilities(); 860 905 int reason2 = 0; … … 2808 2853 ifValuesPass); 2809 2854 #ifdef CLP_USER_DRIVEN 2855 // user can tell which use it is 2856 int status = eventHandler_>event(ClpEventHandler::pivotRow); 2857 if (status >= 0) { 2858 problemStatus_ = 5; 2859 secondaryStatus_ = ClpEventHandler::pivotRow; 2860 break; 2861 } 2810 2862 } else { 2811 2863 int status = eventHandler_>event(ClpEventHandler::pivotRow); … … 2916 2968 userChoiceWasGood(this); 2917 2969 #endif 2918 if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {2970 if (solveType_ >= 2 && (moreSpecialOptions_ & 512) == 0) { 2919 2971 // **** Coding for user interface 2920 2972 // do ray 2921 primalRay(rowArray_[1]); 2973 if (solveType_==2) 2974 primalRay(rowArray_[1]); 2922 2975 // update duals 2923 2976 // as packed need to find pivot row … … 2928 2981 CoinAssert (fabs(alpha_) > 1.0e12); 2929 2982 double multiplier = dualIn_ / alpha_; 2983 #ifndef NDEBUG 2984 rowArray_[0]>checkClear(); 2985 #endif 2930 2986 rowArray_[0]>insert(pivotRow_, multiplier); 2931 2987 factorization_>updateColumnTranspose(rowArray_[2], rowArray_[0]); … … 3147 3203 //printf("costchange on %d from %g to %g\n",sequenceOut_, 3148 3204 // oldCost,cost_[sequenceOut_]); 3149 if (solveType_ !=2)3205 if (solveType_ < 2) 3150 3206 dj_[sequenceOut_] = cost_[sequenceOut_]  oldCost; // normally updated next iteration 3151 3207 solution_[sequenceOut_] = valueOut_; … … 3318 3374 } else if (lower_[iColumn] < 1.0e20 && upper_[iColumn] > 1.0e20) { 3319 3375 setStatus(iColumn, isFree); 3320 break; 3376 if (fabs(dj_[iColumn])>dualTolerance_) 3377 break; 3321 3378 } else { 3322 3379 break;
Note: See TracChangeset
for help on using the changeset viewer.