Ignore:
Timestamp:
Sep 3, 2002 1:43:51 PM (18 years ago)
Author:
forrest
Message:

Synchronizing

File:
1 edited

Legend:

Unmodified
Added
Removed
  • branches/devel-1/ClpSimplexDual.cpp

    r15 r17  
    278278  // for this we need clean basis so it is after factorize
    279279  gutsOfSolution(rowActivityWork_,columnActivityWork_);
     280
     281  numberFake_ =0; // Number of variables at fake bounds
    280282  changeBounds(true,NULL,objectiveChange);
    281283
     
    334336  }
    335337
     338  assert(!numberFake_||problemStatus_); // all bounds should be okay
    336339  // at present we are leaving factorization around
    337340  // maybe we should empty it
     
    346349  return problemStatus_;
    347350}
    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 */
     359int
    349360ClpSimplexDual::whileIterating()
    350361{
    351362  // status stays at -1 while iterating, >=0 finished, -2 to invert
    352363  // status -3 to go to top without an invert
     364  int returnCode = -1;
    353365  while (problemStatus_==-1) {
    354366#ifdef CLP_DEBUG
     
    470482            rowArray_[1]->clear();
    471483            columnArray_[0]->clear();
     484            returnCode=-2;
    472485            break;
    473486          } else {
     
    525538            dualRowPivot_->unrollWeights();
    526539            problemStatus_=-2; // factorize now
     540            returnCode=-2;
    527541            break;
    528542          }
     
    543557            dualRowPivot_->unrollWeights();
    544558            problemStatus_=-2; // factorize now
     559            returnCode=-2;
    545560            break;
    546561          }
     
    552567        if (updateStatus==1) {
    553568          // slight error
    554           if (factorization_->pivots()>5)
     569          if (factorization_->pivots()>5) {
    555570            problemStatus_=-2; // factorize now
     571            returnCode=-3;
     572          }
    556573        } else if (updateStatus==2) {
    557574          // major error
     
    560577          if (factorization_->pivots()) {
    561578            problemStatus_=-2; // factorize now
     579            returnCode=-2;
    562580            break;
    563581          } else {
     
    613631          // maximum iterations or equivalent
    614632          problemStatus_= 3;
     633          returnCode=3;
    615634          break;
    616635        }
     
    630649        rowArray_[0]->clear();
    631650        columnArray_[0]->clear();
     651        returnCode=1;
    632652        break;
    633653      }
     
    660680        }
    661681      }
     682      returnCode=0;
    662683      break;
    663684    }
    664685  }
     686  return returnCode;
    665687}
    666688/* The duals are updated by the given arrays.
     
    838860    // do actual flips
    839861    flipBounds(rowArray,columnArray,theta);
     862    numberFake_ = numberAtFake;
    840863  }
    841864  objectiveChange += changeObj;
     
    17061729    // save dual weights
    17071730    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      }
    17261751    }
    17271752    problemStatus_=-3;
     
    19331958}
    19341959/* While updateDualsInDual sees what effect is of flip
    1935    this does actuall flipping.
     1960   this does actual flipping.
    19361961   If change >0.0 then value in array >0.0 => from lower to upper
    19371962*/
     
    21292154
    21302155}
    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*/
     2162int 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.0e-13);
     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)
     2435int 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_||returnCode||problemStatus_); // 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 */
     2536int
     2537ClpSimplexDual::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::upperFake||bound==ClpSimplex::bothFake)
     2553        numberFake++;
     2554      break;
     2555    case ClpSimplex::atLowerBound:
     2556      if (bound==ClpSimplex::lowerFake||bound==ClpSimplex::bothFake)
     2557        numberFake++;
     2558      break;
     2559    }
     2560  }
     2561  return numberFake;
     2562}
Note: See TracChangeset for help on using the changeset viewer.