Changeset 1761 for trunk


Ignore:
Timestamp:
Jul 6, 2011 12:06:24 PM (8 years ago)
Author:
forrest
Message:

more flexibility in pivoting

Location:
trunk/Clp/src
Files:
9 edited

Legend:

Unmodified
Added
Removed
  • trunk/Clp/src/ClpEventHandler.hpp

    r1665 r1761  
    3434     enum Event {
    3535          endOfIteration = 100, // used to set secondary status
    36           endOfFactorization,
     36          endOfFactorization, // after gutsOfSolution etc
    3737          endOfValuesPass,
    3838          node, // for Cbc
     
    4747          presolveAfterFirstSolve, // ClpSolve presolve after solve
    4848          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)
    5060     };
    5161     /**@name Virtual method that the derived classes should provide.
  • trunk/Clp/src/ClpNonLinearCost.cpp

    r1723 r1761  
    16541654     }
    16551655     if (CLP_METHOD2) {
    1656           abort(); // may never work
     1656       bound_[sequence]=0.0;
     1657       cost2_[sequence]=costValue;
     1658       setInitialStatus(status_[sequence]);
    16571659     }
    16581660
  • trunk/Clp/src/ClpSimplex.cpp

    r1738 r1761  
    80998099     upperIn_ = upper_[sequenceIn_];
    81008100     dualIn_ = dj_[sequenceIn_];
     8101     if (!nonLinearCost_)
     8102       nonLinearCost_ = new ClpNonLinearCost(this);
     8103     
    81018104
    81028105     int returnCode = static_cast<ClpSimplexPrimal *> (this)->pivotResult();
     
    81088111          return -1;
    81098112     }
    8110 }
    8111 
    8112 /* Pivot out a variable and choose an incoing one.  Assumes dual
    8113    feasible - will not go through a reduced cost.
    8114    Returns step length in theta
    8115    Returns ray in ray_ (or NULL if no pivot)
    8116    Return codes as before but -1 means no acceptable pivot
    8117 */
    8118 int
    8119 ClpSimplex::dualPivotResult()
    8120 {
    8121      return static_cast<ClpSimplexDual *> (this)->pivotResult();
    81228113}
    81238114// Factorization frequency
  • trunk/Clp/src/ClpSimplex.hpp

    r1729 r1761  
    340340                       double * valueIncrease, int * sequenceIncrease,
    341341                       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 non-zero 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);
    342360     /** Write the basis in MPS format to the specified file.
    343361         If writeValues true writes values of structurals
     
    458476         feasible - will not go through a reduced cost.
    459477         Returns step length in theta
    460          Returns ray in ray_ (or NULL if no pivot)
    461478         Return codes as before but -1 means no acceptable pivot
    462479     */
    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);
    464485
    465486     /** Common bits of coding for dual and primal.  Return 0 if okay,
     
    839860          return dualIn_;
    840861     }
     862     /// Set reduced cost of last incoming to force error
     863     inline void setDualIn(double value) {
     864          dualIn_ = value;
     865     }
    841866     /// Pivot Row for use by classes e.g. steepestedge
    842867     inline int pivotRow() const {
     
    9821007          valueOut_ = value;
    9831008     }
     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     }
    9841017     /// Set lower of out variable
    9851018     inline void setLowerOut(double value) {
     
    11661199     inline int progressFlag() const {
    11671200          return (progressFlag_ & 3);
     1201     }
     1202     /// For dealing with all issues of cycling etc
     1203     inline ClpSimplexProgress * progress()
     1204     { return &progress_;}
     1205     /// Force re-factorization early value
     1206     inline int forceFactorization() const {
     1207          return forceFactorization_ ;
    11681208     }
    11691209     /// Force re-factorization early
  • trunk/Clp/src/ClpSimplexDual.cpp

    r1729 r1761  
    64436443   feasible - will not go through a reduced cost.
    64446444   Returns step length in theta
    6445    Returns ray in ray_ (or NULL if no pivot)
    64466445   Return codes as before but -1 means no acceptable pivot
    64476446*/
    64486447int
    6449 ClpSimplexDual::pivotResult()
     6448ClpSimplexDual::pivotResultPart1()
    64506449{
    6451      abort();
    6452      return 0;
     6450  // Get good size for pivot
     6451  // Allow first few iterations to take tiny
     6452  double acceptablePivot = 1.0e-1 * 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.0e-8
     6462  //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
     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.0e-8;
     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.0e-6)
     6505    return -1;
     6506  else
     6507    return 0;
    64536508}
    64546509/*
  • trunk/Clp/src/ClpSimplexDual.hpp

    r1665 r1761  
    282282     /** Pivot in a variable and choose an outgoing one.  Assumes dual
    283283         feasible - will not go through a reduced cost.  Returns step length in theta
    284          Returns ray in ray_ (or NULL if no pivot)
    285284         Return codes as before but -1 means no acceptable pivot
    286285     */
    287      int pivotResult();
     286     int pivotResultPart1();
    288287     /** Get next free , -1 if none */
    289288     int nextSuperBasic();
  • trunk/Clp/src/ClpSimplexOther.cpp

    r1738 r1761  
    1616#include "ClpFactorization.hpp"
    1717#include "ClpDualRowDantzig.hpp"
     18#include "ClpNonLinearCost.hpp"
    1819#include "ClpDynamicMatrix.hpp"
    1920#include "CoinPackedMatrix.hpp"
     
    18891890          ClpDualRowPivot * savePivot = dualRowPivot_;
    18901891          dualRowPivot_ = new ClpDualRowDantzig();
     1892          dualRowPivot_->setModel(this);
    18911893
    18921894          if (!returnCode) {
     
    19831985               reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
    19841986               assert (!problemStatus_);
     1987               for (int i=0;i<numberRows_+numberColumns_;i++)
     1988                 setFakeBound(i, noFake);
    19851989               // Now do parametrics
    19861990               handler_->message(CLP_PARAMETRICS_STATS, messages_)
     
    20922096       << line << CoinMessageEol;
    20932097     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*/
     2108int
     2109ClpSimplexOther::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_;
    20942291}
    20952292/* Version of parametrics which reads from file
     
    27692966     }
    27702967}
     2968int
     2969ClpSimplexOther::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>=endingTheta-1.0e-7 )
     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}
    27713075/* Checks if finished.  Updates status */
    27723076void
     
    29143218     // status -3 to go to top without an invert
    29153219     int returnCode = -1;
    2916      double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
    29173220     double lastTheta = startingTheta;
    29183221     double useTheta = startingTheta;
     
    29293232          }
    29303233     }
     3234#if 0
    29313235     // See if objective
    29323236     for (iSequence = 0; iSequence < numberTotal; iSequence++) {
     
    29363240          }
    29373241     }
    2938      assert (type);
     3242#endif
     3243     if (!type) {
     3244       problemStatus_=0;
     3245       return 0;
     3246     }
    29393247     while (problemStatus_ == -1) {
    29403248          double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
     
    29463254          double change = useTheta - lastTheta;
    29473255          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              }
    29633273            }
     3274#if 0
    29643275            cost_[i] += change * changeObjective[i];
     3276#endif
    29653277            assert (solution_[i]>lower_[i]-1.0e-5&&
    29663278                    solution_[i]<upper_[i]+1.0e-5);
     
    29693281          if (pivotType) {
    29703282              problemStatus_ = -2;
     3283              if (!factorization_->pivots()&&pivotRow_<0)
     3284                problemStatus_=2;
    29713285              endingTheta = useTheta;
    29723286              return 4;
     
    29833297               // check accuracy of weights
    29843298               dualRowPivot_->checkAccuracy();
    2985                // Get good size for pivot
    2986                // Allow first few iterations to take tiny
    2987                double acceptablePivot = 1.0e-9;
    2988                if (numberIterations_ > 100)
    2989                     acceptablePivot = 1.0e-8;
    2990                if (factorization_->pivots() > 10 ||
    2991                          (factorization_->pivots() && saveSumDual))
    2992                     acceptablePivot = 1.0e-5; // if we have iterated be more strict
    2993                else if (factorization_->pivots() > 5)
    2994                     acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
    2995                else if (factorization_->pivots())
    2996                     acceptablePivot = 1.0e-8; // relax
    2997                double bestPossiblePivot = 1.0;
    2998                // get sign for finding row of tableau
    2999                // normal iteration
    3000                // create as packed
    3001                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]);
    30073299               // 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();
    30113301               if (sequenceIn_ >= 0) {
    30123302                    // normal iteration
     
    32063496                    objectiveChange -= objectiveValue_;
    32073497                    // outgoing
     3498                    originalBound(sequenceOut_,useTheta,changeLower,changeUpper);
     3499                    lowerOut_=lower_[sequenceOut_];
     3500                    upperOut_=upper_[sequenceOut_];
    32083501                    // set dj to zero unless values pass
    32093502                    if (directionOut_ > 0) {
     
    32543547                    }
    32553548                    // and set bounds correctly
    3256                     reinterpret_cast<ClpSimplexDual *> ( this)->originalBound(sequenceIn_);
     3549                    originalBound(sequenceIn_,useTheta,changeLower,changeUpper);
    32573550                    reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
    32583551                    if (whatNext == 1) {
     
    32763569               } else {
    32773570                    // no incoming column is valid
     3571                 int savePivotRow = pivotRow_;
    32783572                    pivotRow_ = -1;
    32793573#ifdef CLP_DEBUG
     
    32813575                         printf("** no column pivot\n");
    32823576#endif
    3283                     if (factorization_->pivots() < 5) {
     3577                    if (factorization_->pivots() < 5) {
     3578#if 0
    32843579                         // If not in branch and bound etc save ray
    32853580                         if ((specialOptions_&(1024 | 4096)) == 0) {
     
    32903585                              ClpDisjointCopyN(rowArray_[0]->denseVector(), numberRows_, ray_);
    32913586                         }
     3587#endif
    32923588                         // If we have just factorized and infeasibility reasonable say infeas
    32933589                         if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
     
    33033599                                   rowArray_[0]->clear();
    33043600                                   columnArray_[0]->clear();
     3601                                   pivotRow_ = savePivotRow;
     3602                                   int action = eventHandler_->event(ClpEventHandler::noTheta);
     3603                                   if (action==-2)
     3604                                     problemStatus_=-1; // carry on
    33053605                                   returnCode = 1;
    33063606                                   break;
     
    34533753     delete [] primalChange;
    34543754     delete [] dualChange;
    3455      endingTheta = lastTheta;
     3755     endingTheta = lastTheta+theta_;
    34563756     return returnCode;
     3757}
     3758// Finds best possible pivot
     3759double
     3760ClpSimplexOther::bestPivot(bool justColumns)
     3761{
     3762  // Get good size for pivot
     3763  // Allow first few iterations to take tiny
     3764  double acceptablePivot = 1.0e-9;
     3765  if (numberIterations_ > 100)
     3766    acceptablePivot = 1.0e-8;
     3767  if (factorization_->pivots() > 10 ||
     3768      (factorization_->pivots() && sumDualInfeasibilities_))
     3769    acceptablePivot = 1.0e-5; // if we have iterated be more strict
     3770  else if (factorization_->pivots() > 5)
     3771    acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
     3772  else if (factorization_->pivots())
     3773    acceptablePivot = 1.0e-8; // 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;
    34573794}
    34583795// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
     
    35693906          return 0;
    35703907     } else {
     3908          theta_=0.0;
    35713909          return -1;
     3910     }
     3911}
     3912// Restores bound to original bound
     3913void
     3914ClpSimplexOther::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          }
    35723954     }
    35733955}
     
    50065388  //printf("objective value is %g\n", objValue);
    50075389}
     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 non-zero 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*/
     5405int
     5406ClpSimplex::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.0e-7)
     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*/
     5635int
     5636ClpSimplex::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*/
     5644int
     5645ClpSimplex::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.0e-7;
     5821    // if can't trust much and long way from optimal then relax
     5822    if (largestPrimalError_ > 10.0)
     5823      checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
     5824    if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
     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.0e-8 || fabs(alpha_) < 1.0e-8)
     5843          test = 1.0e-1 * fabs(alpha_);
     5844        else
     5845          test = 1.0e-4 * (1.0 + fabs(alpha_));
     5846        if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
     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.0e-60);
     5859    //columnArray_[0]->cleanAndPackSafe(1.0e-60);
     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.0e-5, 1.0e-12 * 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.0e-5)
     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  
    8989         Returns -2 if unable to open file */
    9090     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);
    91105
    92106private:
     
    103117                         const double * changeObjective, ClpDataSave & data,
    104118                         bool canTryQuick);
     119     int parametricsLoop(double startingTheta, double & endingTheta,
     120                         ClpDataSave & data);
    105121     /**  Refactorizes if necessary
    106122          Checks if finished.  Updates status.
     
    131147                   const double * changeLower, const double * changeUpper,
    132148                   const double * changeObjective);
     149     /// Restores bound to original bound
     150     void originalBound(int iSequence, double theta, const double * changeLower,
     151                     const double * changeUpper);
    133152     /**
    134153         Row array has row part of pivot row
  • trunk/Clp/src/ClpSimplexPrimal.cpp

    r1738 r1761  
    460460
    461461               // 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; // re-factorize
     471                 } else if (status>=10) {
     472                   problemStatus_=status-10;
     473                   break;
     474                 } else {
     475                   break;
     476                 }
     477#else
     478                 break;
     479#endif
     480               }
    464481
    465482               // test for maximum iterations
     
    720737               // Force to re-factorize early next time
    721738               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_=status-10;
     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
    725761          }
    726762     }
     
    856892                                            static_cast<double>(numberDualInfeasibilities_ + 10);
    857893#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
    858902          numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
     903#endif
    859904          double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
    860905          int reason2 = 0;
     
    28082853                         ifValuesPass);
    28092854#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               }
    28102862          } else {
    28112863               int status = eventHandler_->event(ClpEventHandler::pivotRow);
     
    29162968               userChoiceWasGood(this);
    29172969 #endif
    2918                if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
     2970               if (solveType_ >= 2 && (moreSpecialOptions_ & 512) == 0) {
    29192971                    // **** Coding for user interface
    29202972                    // do ray
    2921                     primalRay(rowArray_[1]);
     2973                    if (solveType_==2)
     2974                      primalRay(rowArray_[1]);
    29222975                    // update duals
    29232976                    // as packed need to find pivot row
     
    29282981                    CoinAssert (fabs(alpha_) > 1.0e-12);
    29292982                    double multiplier = dualIn_ / alpha_;
     2983#ifndef NDEBUG
     2984                    rowArray_[0]->checkClear();
     2985#endif
    29302986                    rowArray_[0]->insert(pivotRow_, multiplier);
    29312987                    factorization_->updateColumnTranspose(rowArray_[2], rowArray_[0]);
     
    31473203               //printf("costchange on %d from %g to %g\n",sequenceOut_,
    31483204               //       oldCost,cost_[sequenceOut_]);
    3149                if (solveType_ != 2)
     3205               if (solveType_ < 2)
    31503206                    dj_[sequenceOut_] = cost_[sequenceOut_] - oldCost; // normally updated next iteration
    31513207               solution_[sequenceOut_] = valueOut_;
     
    33183374                              } else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
    33193375                                   setStatus(iColumn, isFree);
    3320                                    break;
     3376                                   if (fabs(dj_[iColumn])>dualTolerance_)
     3377                                     break;
    33213378                              } else {
    33223379                                   break;
Note: See TracChangeset for help on using the changeset viewer.