Changeset 17 for branches/devel1/ClpSimplexDual.cpp
 Timestamp:
 Sep 3, 2002 1:43:51 PM (18 years ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel1/ClpSimplexDual.cpp
r15 r17 278 278 // for this we need clean basis so it is after factorize 279 279 gutsOfSolution(rowActivityWork_,columnActivityWork_); 280 281 numberFake_ =0; // Number of variables at fake bounds 280 282 changeBounds(true,NULL,objectiveChange); 281 283 … … 334 336 } 335 337 338 assert(!numberFake_problemStatus_); // all bounds should be okay 336 339 // at present we are leaving factorization around 337 340 // maybe we should empty it … … 346 349 return problemStatus_; 347 350 } 348 void 351 /* Reasons to come out: 352 1 iterations etc 353 2 inaccuracy 354 3 slight inaccuracy (and done iterations) 355 +0 looks optimal (might be unbounded  but we will investigate) 356 +1 looks infeasible 357 +3 max iterations 358 */ 359 int 349 360 ClpSimplexDual::whileIterating() 350 361 { 351 362 // status stays at 1 while iterating, >=0 finished, 2 to invert 352 363 // status 3 to go to top without an invert 364 int returnCode = 1; 353 365 while (problemStatus_==1) { 354 366 #ifdef CLP_DEBUG … … 470 482 rowArray_[1]>clear(); 471 483 columnArray_[0]>clear(); 484 returnCode=2; 472 485 break; 473 486 } else { … … 525 538 dualRowPivot_>unrollWeights(); 526 539 problemStatus_=2; // factorize now 540 returnCode=2; 527 541 break; 528 542 } … … 543 557 dualRowPivot_>unrollWeights(); 544 558 problemStatus_=2; // factorize now 559 returnCode=2; 545 560 break; 546 561 } … … 552 567 if (updateStatus==1) { 553 568 // slight error 554 if (factorization_>pivots()>5) 569 if (factorization_>pivots()>5) { 555 570 problemStatus_=2; // factorize now 571 returnCode=3; 572 } 556 573 } else if (updateStatus==2) { 557 574 // major error … … 560 577 if (factorization_>pivots()) { 561 578 problemStatus_=2; // factorize now 579 returnCode=2; 562 580 break; 563 581 } else { … … 613 631 // maximum iterations or equivalent 614 632 problemStatus_= 3; 633 returnCode=3; 615 634 break; 616 635 } … … 630 649 rowArray_[0]>clear(); 631 650 columnArray_[0]>clear(); 651 returnCode=1; 632 652 break; 633 653 } … … 660 680 } 661 681 } 682 returnCode=0; 662 683 break; 663 684 } 664 685 } 686 return returnCode; 665 687 } 666 688 /* The duals are updated by the given arrays. … … 838 860 // do actual flips 839 861 flipBounds(rowArray,columnArray,theta); 862 numberFake_ = numberAtFake; 840 863 } 841 864 objectiveChange += changeObj; … … 1706 1729 // save dual weights 1707 1730 dualRowPivot_>saveWeights(this,1); 1708 // is factorization okay? 1709 if (internalFactorize(1)) { 1710 // no  restore previous basis 1711 assert (type==1); 1712 changeMade_++; // say something changed 1713 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char)); 1714 memcpy(rowActivityWork_,savedSolution_+numberColumns_ , 1715 numberRows_*sizeof(double)); 1716 memcpy(columnActivityWork_,savedSolution_ , 1717 numberColumns_*sizeof(double)); 1718 // get correct bounds on all variables 1719 double dummyChangeCost=0.0; 1720 changeBounds(true,rowArray_[2],dummyChangeCost); 1721 // throw away change 1722 rowArray_[2]>clear(); 1723 forceFactorization_=1; // a bit drastic but .. 1724 type = 2; 1725 assert (internalFactorize(1)==0); 1731 if (type) { 1732 // is factorization okay? 1733 if (internalFactorize(1)) { 1734 // no  restore previous basis 1735 assert (type==1); 1736 changeMade_++; // say something changed 1737 memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char)); 1738 memcpy(rowActivityWork_,savedSolution_+numberColumns_ , 1739 numberRows_*sizeof(double)); 1740 memcpy(columnActivityWork_,savedSolution_ , 1741 numberColumns_*sizeof(double)); 1742 // get correct bounds on all variables 1743 double dummyChangeCost=0.0; 1744 changeBounds(true,rowArray_[2],dummyChangeCost); 1745 // throw away change 1746 rowArray_[2]>clear(); 1747 forceFactorization_=1; // a bit drastic but .. 1748 type = 2; 1749 assert (internalFactorize(1)==0); 1750 } 1726 1751 } 1727 1752 problemStatus_=3; … … 1933 1958 } 1934 1959 /* While updateDualsInDual sees what effect is of flip 1935 this does actual lflipping.1960 this does actual flipping. 1936 1961 If change >0.0 then value in array >0.0 => from lower to upper 1937 1962 */ … … 2129 2154 2130 2155 } 2131 2156 /* For strong branching. On input lower and upper are new bounds 2157 while on output they are change in objective function values 2158 (>1.0e50 infeasible). 2159 Return code is 0 if nothing interesting, 1 if infeasible both 2160 ways and +1 if infeasible one way (check values to see which one(s)) 2161 */ 2162 int ClpSimplexDual::strongBranching(int numberVariables,const int * variables, 2163 double * newLower, double * newUpper, 2164 bool stopOnFirstInfeasible, 2165 bool alwaysFinish) 2166 { 2167 int i; 2168 int returnCode=0; 2169 double saveObjective = objectiveValue_; 2170 #if 1 2171 algorithm_ = 1; 2172 2173 //scaling(false); 2174 2175 // put in standard form (and make row copy) 2176 // create modifiable copies of model rim and do optional scaling 2177 createRim(7+8+16,true); 2178 2179 // change newLower and newUpper if scaled 2180 2181 // Do initial factorization 2182 // and set certain stuff 2183 // We can either set increasing rows so ...IsBasic gives pivot row 2184 // or we can just increment iBasic one by one 2185 // for now let ...iBasic give pivot row 2186 factorization_>increasingRows(2); 2187 // row activities have negative sign 2188 factorization_>slackValue(1.0); 2189 factorization_>zeroTolerance(1.0e13); 2190 // save if sparse factorization wanted 2191 int saveSparse = factorization_>sparseThreshold(); 2192 2193 int factorizationStatus = internalFactorize(0); 2194 if (factorizationStatus<0) 2195 return 1; // some error 2196 else if (factorizationStatus) 2197 handler_>message(CLP_SINGULARITIES,messages_) 2198 <<factorizationStatus 2199 <<OsiMessageEol; 2200 if (saveSparse) { 2201 // use default at present 2202 factorization_>sparseThreshold(0); 2203 factorization_>goSparse(); 2204 } 2205 2206 // save stuff 2207 ClpFactorization saveFactorization(*factorization_); 2208 // save basis and solution 2209 double * saveSolution = new double[numberRows_+numberColumns_]; 2210 memcpy(saveSolution,solution_, 2211 (numberRows_+numberColumns_)*sizeof(double)); 2212 unsigned char * saveStatus = 2213 new unsigned char [numberRows_+numberColumns_]; 2214 memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char)); 2215 // save bounds as createRim makes clean copies 2216 double * saveLower = new double[numberRows_+numberColumns_]; 2217 memcpy(saveLower,lower_, 2218 (numberRows_+numberColumns_)*sizeof(double)); 2219 double * saveUpper = new double[numberRows_+numberColumns_]; 2220 memcpy(saveUpper,upper_, 2221 (numberRows_+numberColumns_)*sizeof(double)); 2222 int * savePivot = new int [numberRows_]; 2223 memcpy(savePivot, pivotVariable_, numberRows_*sizeof(int)); 2224 // need to save/restore weights. 2225 2226 for (i=0;i<numberVariables;i++) { 2227 int iColumn = variables[i]; 2228 double objectiveChange; 2229 double saveBound; 2230 2231 // try down 2232 2233 saveBound = columnUpper_[iColumn]; 2234 // external view  in case really getting optimal 2235 columnUpper_[iColumn] = newUpper[i]; 2236 if (scalingFlag_<=0) 2237 upper_[iColumn] = newUpper[i]; 2238 else 2239 upper_[iColumn] = newUpper[i]/columnScale_[iColumn]; // scale 2240 // Start of fast iterations 2241 int status = fastDual(alwaysFinish); 2242 2243 // restore 2244 memcpy(solution_,saveSolution, 2245 (numberRows_+numberColumns_)*sizeof(double)); 2246 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char)); 2247 memcpy(lower_,saveLower, 2248 (numberRows_+numberColumns_)*sizeof(double)); 2249 memcpy(upper_,saveUpper, 2250 (numberRows_+numberColumns_)*sizeof(double)); 2251 columnUpper_[iColumn] = saveBound; 2252 memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int)); 2253 delete factorization_; 2254 factorization_ = new ClpFactorization(saveFactorization); 2255 2256 if (status(problemStatus_==0&&!isDualObjectiveLimitReached())) { 2257 objectiveChange = objectiveValue_saveObjective; 2258 } else { 2259 objectiveChange = 1.0e100; 2260 } 2261 newUpper[i]=objectiveChange; 2262 #ifdef CLP_DEBUG 2263 printf("down on %d costs %g\n",iColumn,objectiveChange); 2264 #endif 2265 2266 // try up 2267 2268 saveBound = columnLower_[iColumn]; 2269 // external view  in case really getting optimal 2270 columnLower_[iColumn] = newLower[i]; 2271 if (scalingFlag_<=0) 2272 lower_[iColumn] = newLower[i]; 2273 else 2274 lower_[iColumn] = newLower[i]/columnScale_[iColumn]; // scale 2275 // Start of fast iterations 2276 status = fastDual(alwaysFinish); 2277 2278 // restore 2279 memcpy(solution_,saveSolution, 2280 (numberRows_+numberColumns_)*sizeof(double)); 2281 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char)); 2282 memcpy(lower_,saveLower, 2283 (numberRows_+numberColumns_)*sizeof(double)); 2284 memcpy(upper_,saveUpper, 2285 (numberRows_+numberColumns_)*sizeof(double)); 2286 columnLower_[iColumn] = saveBound; 2287 memcpy(pivotVariable_, savePivot, numberRows_*sizeof(int)); 2288 delete factorization_; 2289 factorization_ = new ClpFactorization(saveFactorization); 2290 2291 if (status(problemStatus_==0&&!isDualObjectiveLimitReached())) { 2292 objectiveChange = objectiveValue_saveObjective; 2293 } else { 2294 objectiveChange = 1.0e100; 2295 } 2296 newLower[i]=objectiveChange; 2297 #ifdef CLP_DEBUG 2298 printf("up on %d costs %g\n",iColumn,objectiveChange); 2299 #endif 2300 2301 /* Possibilities are: 2302 Both sides feasible  store 2303 Neither side feasible  set objective high and exit 2304 One side feasible  change bounds and resolve 2305 */ 2306 if (newUpper[i]<1.0e100) { 2307 if(newLower[i]<1.0e100) { 2308 // feasible  no action 2309 } else { 2310 // up feasible, down infeasible 2311 returnCode=1; 2312 if (stopOnFirstInfeasible) 2313 break; 2314 } 2315 } else { 2316 if(newLower[i]<1.0e100) { 2317 // down feasible, up infeasible 2318 returnCode=1; 2319 if (stopOnFirstInfeasible) 2320 break; 2321 } else { 2322 // neither side feasible 2323 returnCode=1; 2324 break; 2325 } 2326 } 2327 } 2328 delete [] saveSolution; 2329 delete [] saveLower; 2330 delete [] saveUpper; 2331 delete [] saveStatus; 2332 delete [] savePivot; 2333 2334 // at present we are leaving factorization around 2335 // maybe we should empty it 2336 deleteRim(); 2337 factorization_>sparseThreshold(saveSparse); 2338 #else 2339 // save basis and solution 2340 double * rowSolution = new double[numberRows_]; 2341 memcpy(rowSolution,rowActivity_,numberRows_*sizeof(double)); 2342 double * columnSolution = new double[numberColumns_]; 2343 memcpy(columnSolution,columnActivity_,numberColumns_*sizeof(double)); 2344 unsigned char * saveStatus = 2345 new unsigned char [numberRows_+numberColumns_]; 2346 if (!status_) { 2347 // odd but anyway setup all slack basis 2348 status_ = new unsigned char [numberColumns_+numberRows_]; 2349 memset(status_,0,(numberColumns_+numberRows_)*sizeof(char)); 2350 } 2351 memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char)); 2352 for (i=0;i<numberVariables;i++) { 2353 int iColumn = variables[i]; 2354 double objectiveChange; 2355 2356 // try down 2357 2358 double saveUpper = columnUpper_[iColumn]; 2359 columnUpper_[iColumn] = newUpper[i]; 2360 dual(); 2361 2362 // restore 2363 columnUpper_[iColumn] = saveUpper; 2364 memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double)); 2365 memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double)); 2366 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char)); 2367 2368 if (problemStatus_==0&&!isDualObjectiveLimitReached()) { 2369 objectiveChange = objectiveValue_saveObjective; 2370 } else { 2371 objectiveChange = 1.0e100; 2372 } 2373 newUpper[i]=objectiveChange; 2374 #ifdef CLP_DEBUG 2375 printf("down on %d costs %g\n",iColumn,objectiveChange); 2376 #endif 2377 2378 // try up 2379 2380 double saveLower = columnLower_[iColumn]; 2381 columnLower_[iColumn] = newLower[i]; 2382 dual(); 2383 2384 // restore 2385 columnLower_[iColumn] = saveLower; 2386 memcpy(rowActivity_,rowSolution,numberRows_*sizeof(double)); 2387 memcpy(columnActivity_,columnSolution,numberColumns_*sizeof(double)); 2388 memcpy(status_,saveStatus,(numberColumns_+numberRows_)*sizeof(char)); 2389 2390 if (problemStatus_==0&&!isDualObjectiveLimitReached()) { 2391 objectiveChange = objectiveValue_saveObjective; 2392 } else { 2393 objectiveChange = 1.0e100; 2394 } 2395 newLower[i]=objectiveChange; 2396 #ifdef CLP_DEBUG 2397 printf("up on %d costs %g\n",iColumn,objectiveChange); 2398 #endif 2399 2400 /* Possibilities are: 2401 Both sides feasible  store 2402 Neither side feasible  set objective high and exit 2403 One side feasible  change bounds and resolve 2404 */ 2405 if (newUpper[i]<1.0e100) { 2406 if(newLower[i]<1.0e100) { 2407 // feasible  no action 2408 } else { 2409 // up feasible, down infeasible 2410 returnCode=1; 2411 if (stopOnFirstInfeasible) 2412 break; 2413 } 2414 } else { 2415 if(newLower[i]<1.0e100) { 2416 // down feasible, up infeasible 2417 returnCode=1; 2418 if (stopOnFirstInfeasible) 2419 break; 2420 } else { 2421 // neither side feasible 2422 returnCode=1; 2423 break; 2424 } 2425 } 2426 } 2427 delete [] rowSolution; 2428 delete [] columnSolution; 2429 delete [] saveStatus; 2430 #endif 2431 objectiveValue_ = saveObjective; 2432 return returnCode; 2433 } 2434 // treat no pivot as finished (unless interesting) 2435 int ClpSimplexDual::fastDual(bool alwaysFinish) 2436 { 2437 algorithm_ = 1; 2438 dualTolerance_=dblParam_[OsiDualTolerance]; 2439 primalTolerance_=dblParam_[OsiPrimalTolerance]; 2440 2441 // save dual bound 2442 double saveDualBound = dualBound_; 2443 2444 double objectiveChange; 2445 // for dual we will change bounds using dualBound_ 2446 // for this we need clean basis so it is after factorize 2447 gutsOfSolution(rowActivityWork_,columnActivityWork_); 2448 numberFake_ =0; // Number of variables at fake bounds 2449 changeBounds(true,NULL,objectiveChange); 2450 2451 problemStatus_ = 1; 2452 numberIterations_=0; 2453 2454 int lastCleaned=0; // last time objective or bounds cleaned up 2455 2456 // number of times we have declared optimality 2457 numberTimesOptimal_=0; 2458 2459 // This says whether to restore things etc 2460 int factorType=0; 2461 /* 2462 Status of problem: 2463 0  optimal 2464 1  infeasible 2465 2  unbounded 2466 1  iterating 2467 2  factorization wanted 2468 3  redo checking without factorization 2469 4  looks infeasible 2470 2471 BUT also from whileIterating return code is: 2472 2473 1 iterations etc 2474 2 inaccuracy 2475 3 slight inaccuracy (and done iterations) 2476 +0 looks optimal (might be unbounded  but we will investigate) 2477 +1 looks infeasible 2478 +3 max iterations 2479 2480 */ 2481 2482 int returnCode = 0; 2483 2484 while (problemStatus_<0) { 2485 int iRow,iColumn; 2486 // clear 2487 for (iRow=0;iRow<4;iRow++) { 2488 rowArray_[iRow]>clear(); 2489 } 2490 2491 for (iColumn=0;iColumn<2;iColumn++) { 2492 columnArray_[iColumn]>clear(); 2493 } 2494 2495 // give matrix (and model costs and bounds a chance to be 2496 // refreshed (normally null) 2497 matrix_>refresh(this); 2498 // may factorize, checks if problem finished 2499 // should be able to speed this up on first time 2500 statusOfProblemInDual(lastCleaned,factorType); 2501 2502 // Say good factorization 2503 factorType=1; 2504 2505 // Do iterations 2506 if (problemStatus_<0) { 2507 #if 1 2508 returnCode = whileIterating(); 2509 if (!alwaysFinish&&returnCode<1) { 2510 double limit = 0.0; 2511 getDblParam(OsiDualObjectiveLimit, limit); 2512 if(objectiveValue()*optimizationDirection_<limit 2513 numberAtFakeBound()) { 2514 // can't say anything interesting  might as well return 2515 #ifdef CLP_DEBUG 2516 printf("returning from fastDual after %d iterations with code %d\n", 2517 numberIterations_,returnCode); 2518 #endif 2519 returnCode=1; 2520 break; 2521 } 2522 } 2523 returnCode=0; 2524 #else 2525 whileIterating(); 2526 #endif 2527 } 2528 } 2529 2530 assert(!numberFake_returnCodeproblemStatus_); // all bounds should be okay 2531 dualBound_ = saveDualBound; 2532 return returnCode; 2533 } 2534 /* Checks number of variables at fake bounds. This is used by fastDual 2535 so can exit gracefully before end */ 2536 int 2537 ClpSimplexDual::numberAtFakeBound() 2538 { 2539 int iSequence; 2540 int numberFake=0; 2541 2542 for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) { 2543 FakeBound bound = getFakeBound(iSequence); 2544 switch(getStatus(iSequence)) { 2545 2546 case ClpSimplex::basic: 2547 break; 2548 case ClpSimplex::isFree: 2549 case ClpSimplex::superBasic: 2550 break; 2551 case ClpSimplex::atUpperBound: 2552 if (bound==ClpSimplex::upperFakebound==ClpSimplex::bothFake) 2553 numberFake++; 2554 break; 2555 case ClpSimplex::atLowerBound: 2556 if (bound==ClpSimplex::lowerFakebound==ClpSimplex::bothFake) 2557 numberFake++; 2558 break; 2559 } 2560 } 2561 return numberFake; 2562 }
Note: See TracChangeset
for help on using the changeset viewer.