Changeset 1785

Timestamp:

Aug 25, 2011 6:17:49 AM (9 years ago)

Author:

forrest

Message:

patches to try and make parametrics faster

Location:

trunk/Clp/src

Files:

• ClpSimplex.cpp (modified) (2 diffs)
• ClpSimplexDual.cpp (modified) (8 diffs)
• ClpSimplexOther.cpp (modified) (29 diffs)
• ClpSimplexOther.hpp (modified) (5 diffs)

trunk/Clp/src/ClpSimplex.cpp

@@ -669 +669 @@
           }
      }
+#if CAN_HAVE_ZERO_OBJ>1
+     if ((specialOptions_&2097152)==0) {
+#endif
      computeDuals(givenDuals);
      if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
@@ -692 +695 @@
      checkBothSolutions();
      objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
+#if CAN_HAVE_ZERO_OBJ>1
+     } else {
+       checkPrimalSolution( rowActivityWork_, columnActivityWork_);
+#ifndef COIN_REUSE_RANDOM
+       memset(dj_,0,(numberRows_+numberColumns_)*sizeof(double));
+#else
+       for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
+         double value;
+         switch (getStatus(iSequence)) {
+         case atLowerBound:
+           value=1.0e-9*(1.0+CoinDrand48());
+           break;
+         case atUpperBound:
+           value=-1.0e-9*(1.0+CoinDrand48());
+           break;
+         default:
+           value=0.0;
+           break;
+         }
+       }
+#endif
+       objectiveValue_=0.0;
+     }
+#endif
      if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 ||
                                       largestDualError_ > 1.0e-2))

trunk/Clp/src/ClpSimplexDual.cpp

@@ -1223 +1223 @@
                     bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
                                                    columnArray_[1], acceptablePivot, dubiousWeights);
+#if CAN_HAVE_ZERO_OBJ>1
+                    if ((specialOptions_&2097152)!=0) 
+                      theta_=0.0;
+#endif
                } else {
                     // Make sure direction plausible
@@ -1399 +1403 @@
                     //columnArray_[0]->cleanAndPackSafe(1.0e-60);
                     if (candidate == -1) {
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    if ((specialOptions_&2097152)==0) {
+#endif
                          // make sure incoming doesn't count
                          Status saveStatus = getStatus(sequenceIn_);
@@ -1406 +1413 @@
                                                       objectiveChange, false);
                          setStatus(sequenceIn_, saveStatus);
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    } else {
+                      rowArray_[0]->clear();
+                      rowArray_[2]->clear();
+                      columnArray_[0]->clear();
+                    }
+#endif
                     } else {
                          updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
@@ -3778 +3792 @@
                               if (value > dualTolerance_) {
                                    thisIncrease = true;
+#if CLP_CAN_HAVE_ZERO_OBJ<2
 #define MODIFYCOST 2
+#endif
 #if MODIFYCOST
                                    // modify cost to hit new tolerance
@@ -4666 +4682 @@
                          if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
                                    (numberChanged_ || (specialOptions_ & 4096) == 0)) {
+#if CLP_CAN_HAVE_ZERO_OBJ
                            if ((specialOptions_&2097152)==0) {
+#endif
                               doOriginalTolerance = 2;
                               numberTimesOptimal_++;
@@ -4681 +4699 @@
                               }
                               cleanDuals = 2; // If nothing changed optimal else primal
+#if CLP_CAN_HAVE_ZERO_OBJ
                            } else {
                              // no cost - skip checks
                              problemStatus_=0;
                            }
+#endif
                          } else {
                               problemStatus_ = 0; // optimal
@@ -4939 +4959 @@
                     if (givenDuals)
                          dualTolerance_ = 1.0e50;
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    if ((specialOptions_&2097152)==0) {
+#endif
                     updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
                                       0.0, objectiveChange, true);
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    } else {
+                      rowArray_[0]->clear();
+                      rowArray_[1]->clear();
+                      columnArray_[0]->clear();
+                    }
+#endif
                     dualTolerance_ = saveTolerance;
 #if 0
@@ -4971 +5001 @@
                     if (givenDuals)
                          dualTolerance_ = 1.0e50;
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    if ((specialOptions_&2097152)==0) {
+#endif
                     updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
                                       0.0, objectiveChange, true);
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+                    } else {
+                      rowArray_[0]->clear();
+                      rowArray_[1]->clear();
+                      columnArray_[0]->clear();
+                    }
+#endif
                     dualTolerance_ = saveTolerance;
                     if (!numberIterations_ && sumPrimalInfeasibilities_ >

trunk/Clp/src/ClpSimplexOther.cpp

@@ -942 +942 @@
 // Restores solution from dualized problem
 int
-ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
+ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem,
+                                 bool checkAccuracy)
 {
      int returnCode = 0;;
@@ -1187 +1188 @@
      // Below will go to ..DEBUG later
 #if 1 //ndef NDEBUG
-     // Check if correct
-     double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
-     double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
-     double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
-     double * dual = CoinCopyOfArray(dual_, numberRows_);
-     this->dual(); //primal();
-     CoinRelFltEq eq(1.0e-5);
-     for (iRow = 0; iRow < numberRows_; iRow++) {
-          assert(eq(dual[iRow], dual_[iRow]));
+     if (checkAccuracy) {
+       // Check if correct
+       double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
+       double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
+       double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
+       double * dual = CoinCopyOfArray(dual_, numberRows_);
+       this->dual(); //primal();
+       CoinRelFltEq eq(1.0e-5);
+       for (iRow = 0; iRow < numberRows_; iRow++) {
+         assert(eq(dual[iRow], dual_[iRow]));
+       }
+       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
+         assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
+       }
+       for (iRow = 0; iRow < numberRows_; iRow++) {
+         assert(eq(rowActivity[iRow], rowActivity_[iRow]));
+       }
+       for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
+         assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
+       }
+       delete [] columnActivity;
+       delete [] rowActivity;
+       delete [] reducedCost;
+       delete [] dual;
      }
-     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
-          assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
-     }
-     for (iRow = 0; iRow < numberRows_; iRow++) {
-          assert(eq(rowActivity[iRow], rowActivity_[iRow]));
-     }
-     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
-          assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
-     }
-     delete [] columnActivity;
-     delete [] rowActivity;
-     delete [] reducedCost;
-     delete [] dual;
 #endif
      return returnCode;
@@ -1992 +1995 @@
                while (!returnCode) {
                     //assert (reportIncrement);
-                    returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
+                 parametricsData paramData;
+                 paramData.startingTheta=startingTheta;
+                 paramData.endingTheta=endingTheta;
+                 paramData.lowerChange = chgLower;
+                 paramData.upperChange = chgUpper;
+                    returnCode = parametricsLoop(paramData, reportIncrement,
                                                  chgLower, chgUpper, chgObjective, data,
                                                  canTryQuick);
+                 startingTheta=paramData.startingTheta;
+                 endingTheta=paramData.endingTheta;
                     if (!returnCode) {
                          //double change = endingTheta-startingTheta;
@@ -2656 +2666 @@
 }
 int
-ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
+ClpSimplexOther::parametricsLoop(parametricsData & paramData,double reportIncrement,
                                  const double * lowerChange, const double * upperChange,
                                  const double * changeObjective, ClpDataSave & data,
                                  bool canTryQuick)
 {
+  double startingTheta = paramData.startingTheta;
+  double & endingTheta = paramData.endingTheta;
      // stuff is already at starting
      // For this crude version just try and go to end
@@ -2757 +2769 @@
                reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
           } else {
-               whileIterating(startingTheta,  endingTheta, reportIncrement,
-                              lowerChange, upperChange,
+               whileIterating(paramData, reportIncrement,
                               changeObjective);
                startingTheta = endingTheta;
@@ -2804 +2815 @@
   assert (ratio==1||ratio==2);
   // allow for unscaled - even if not needed
-  int lengthArrays = 4*numberTotal+(numberTotal+2)*ratio;
+  int lengthArrays = 4*numberTotal+(2*numberTotal+2)*ratio;
   /*
     Save information and modify
@@ -2819 +2830 @@
   char * markDone = reinterpret_cast<char *>(lowerList+numberTotal);
   memset(markDone,0,numberTotal);
+  int * backwardBasic = lowerList+2*numberTotal;
+  parametricsData paramData;
+  paramData.lowerChange = lowerChange;
+  paramData.lowerList=lowerList;
+  paramData.upperChange = upperChange;
+  paramData.upperList=upperList;
+  paramData.markDone=markDone;
+  paramData.backwardBasic=backwardBasic;
   // Find theta when bounds will cross over and create arrays
   memset(lowerChange, 0, numberTotal * sizeof(double));
@@ -3036 +3055 @@
         bool canSkipFactorization=true;
         while (!returnCode) {
-          returnCode = parametricsLoop(startingTheta, endingTheta,
+                 paramData.startingTheta=startingTheta;
+                 paramData.endingTheta=endingTheta;
+                 returnCode = parametricsLoop(paramData,
                                        data,canSkipFactorization);
+                 startingTheta=paramData.startingTheta;
+                 endingTheta=paramData.endingTheta;
           canSkipFactorization=false;
           if (!returnCode) {
@@ -3064 +3087 @@
       ray_ = new double [numberColumns_];
     }
-    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(1);
     if (swapped&&lower_) {
       double * temp=saveLower;
@@ -3073 +3095 @@
       upper_=temp;
     }
+    reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
   }    
   if (!scalingFlag_) {
@@ -3117 +3140 @@
 }
 int
-ClpSimplexOther::parametricsLoop(double & startingTheta, double & endingTheta,
+ClpSimplexOther::parametricsLoop(parametricsData & paramData,
                                  ClpDataSave & data,bool canSkipFactorization)
 {
+  double & startingTheta = paramData.startingTheta;
+  double & endingTheta = paramData.endingTheta;
   int numberTotal = numberRows_+numberColumns_;
-  const double * lowerChange = lower_+numberTotal;
-  const double * upperChange = upper_+numberTotal;
   // stuff is already at starting
-  int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
-  int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
+  const int * lowerList = paramData.lowerList;
+  const int * upperList = paramData.upperList;
   problemStatus_ = -1;
   //double saveEndingTheta=endingTheta;
@@ -3259 +3282 @@
     // Do iterations
     problemStatus_=-1;
-    whileIterating(startingTheta,  endingTheta, 0.0,
-                   lowerChange, upperChange,
+    whileIterating(paramData, 0.0,
                    NULL);
     //startingTheta = endingTheta;
@@ -3417 +3439 @@
 */
 int
-ClpSimplexOther::whileIterating(double & startingTheta, double & endingTheta, double /*reportIncrement*/,
-                                const double * lowerChange, const double * upperChange,
+ClpSimplexOther::whileIterating(parametricsData & paramData, double /*reportIncrement*/,
                                 const double * /*changeObjective*/)
 {
+  double & startingTheta = paramData.startingTheta;
+  double & endingTheta = paramData.endingTheta;
+  const double * lowerChange = paramData.lowerChange;
+  const double * upperChange = paramData.upperChange;
   int numberTotal = numberColumns_ + numberRows_;
-  const int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
-  const int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
+  const int * lowerList = paramData.lowerList;
+  const int * upperList = paramData.upperList;
+  // do basic pointers
+  int * backwardBasic = paramData.backwardBasic;
+  for (int i=0;i<numberTotal;i++)
+    backwardBasic[i]=-1;
+  for (int i=0;i<numberRows_;i++) {
+    int iPivot=pivotVariable_[i];
+    backwardBasic[iPivot]=i;
+  }
      {
           int i;
@@ -3446 +3479 @@
           double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
           // Get theta for bounds - we know can't crossover
-          int pivotType = nextTheta(1, increaseTheta, 
-                                    lowerChange, upperChange, NULL);
+          int pivotType = nextTheta(1, increaseTheta, paramData,
+                                     NULL);
           useTheta += theta_;
           double change = useTheta - lastTheta;
@@ -3591 +3624 @@
                     //rowArray_[0]->cleanAndPackSafe(1.0e-60);
                     //columnArray_[0]->cleanAndPackSafe(1.0e-60);
+#if CLP_CAN_HAVE_ZERO_OBJ
                     if ((specialOptions_&2097152)==0) {
+#endif
                       nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
                                                                                                rowArray_[2], theta_,
                                                                                                objectiveChange, false);
                       assert (!nswapped);
+#if CLP_CAN_HAVE_ZERO_OBJ
                     } else {
                       rowArray_[0]->clear();
@@ -3601 +3637 @@
                       columnArray_[0]->clear();
                     }
+#endif
                     // which will change basic solution
                     if (nswapped) {
@@ -3710 +3747 @@
                     }
                     // update primal solution
+#if CLP_CAN_HAVE_ZERO_OBJ
+                    if ((specialOptions_&2097152)!=0) 
+                      theta_=0.0;
+#endif
                     if (theta_ < 0.0) {
 #ifdef CLP_DEBUG
@@ -3737 +3778 @@
                     }
                     objectiveChange = 0.0;
+#if CLP_CAN_HAVE_ZERO_OBJ
                     if ((specialOptions_&2097152)==0) {
+#endif
                       for (int i=0;i<numberTotal;i++)
                         objectiveChange += solution_[i]*cost_[i];
                       objectiveChange -= objectiveValue_;
+#if CLP_CAN_HAVE_ZERO_OBJ
                     }
+#endif
                     // outgoing
                     originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
@@ -3750 +3795 @@
                          valueOut_ = lowerOut_;
                          dj_[sequenceOut_] = theta_;
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+#ifdef COIN_REUSE_RANDOM
+                         if ((specialOptions_&2097152)!=0) {
+                           dj_[sequenceOut_] = 1.0e-9*(1.0+CoinDrand48());;
+                         }
+#endif
+#endif
                     } else {
                          valueOut_ = upperOut_;
                          dj_[sequenceOut_] = -theta_;
+#if CLP_CAN_HAVE_ZERO_OBJ>1
+#ifdef COIN_REUSE_RANDOM
+                         if ((specialOptions_&2097152)!=0) {
+                           dj_[sequenceOut_] = -1.0e-9*(1.0+CoinDrand48());;
+                         }
+#endif
+#endif
                     }
                     solution_[sequenceOut_] = valueOut_;
                     int whatNext = housekeeping(objectiveChange);
+                    assert (backwardBasic[sequenceOut_]==pivotRow_);
+                    backwardBasic[sequenceOut_]=-1;
+                    backwardBasic[sequenceIn_]=pivotRow_;
                     {
                       char in[200],out[200];
@@ -3828 +3890 @@
                    theta_ = 0.0;
                    //adjust [4] from handler - but
-                   rowArray_[4]->clear(); // temp
+                   //rowArray_[4]->clear(); // temp
                    if (action>=0&&action<10)
                      problemStatus_=-1; // carry on
@@ -4063 +4125 @@
 // Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
 int
-ClpSimplexOther::nextTheta(int type, double maxTheta, 
-                           const double * lowerChange, const double * upperChange,
+ClpSimplexOther::nextTheta(int type, double maxTheta, parametricsData & paramData,
                            const double * /*changeObjective*/)
 {
-  int numberTotal = numberColumns_ + numberRows_;
-  const int * lowerList = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
-  const int * upperList = (reinterpret_cast<int *>(upper_+4*numberTotal))+2;
+  const double * lowerChange = paramData.lowerChange;
+  const double * upperChange = paramData.upperChange;
+  const int * lowerList = paramData.lowerList;
+  const int * upperList = paramData.upperList;
   int iSequence;
   theta_ = maxTheta;
@@ -4081 +4143 @@
   const CoinBigIndex * columnStart = matrix_->getVectorStarts();
   const double * elementByColumn = matrix_->getElements();
+#if 0
+  double tempArray[5000];
+  bool checkIt=false;
+  if (factorization_->pivots()&&!needFullUpdate&&sequenceIn_<0) {
+    memcpy(tempArray,array,numberRows_*sizeof(double));
+    checkIt=true;
+    needFullUpdate=true;
+  }
+#endif
   if (!factorization_->pivots()||needFullUpdate) {
     rowArray_[4]->clear();
@@ -4187 +4258 @@
     // ftran it
     factorization_->updateColumn(rowArray_[0], rowArray_[4]);
-  } else {
-    assert (sequenceIn_>=0);
+#if 0
+    if (checkIt) {
+      for (int i=0;i<numberRows_;i++) {
+        assert (fabs(tempArray[i]-array[i])<1.0e-8);
+      }
+    }
+#endif
+  } else if (sequenceIn_>=0) {
+    //assert (sequenceIn_>=0);
     assert (sequenceOut_>=0);
     assert (sequenceIn_!=sequenceOut_);
@@ -4260 +4338 @@
   const int * index = rowArray_[4]->getIndices();
   int number = rowArray_[4]->getNumElements();
-  int * lowerList2 = (reinterpret_cast<int *>(lower_+4*numberTotal))+2;
-  char * markDone = reinterpret_cast<char *>(lowerList2+numberTotal);
-#if 0
-  for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
-    //int iPivot = index[iRow];
-    iSequence = pivotVariable_[iPivot];
-    // solution value will be sol - theta*alpha
-    // bounds will be bounds + change *theta
-    double currentSolution = solution_[iSequence];
-    double currentLower = lower_[iSequence];
-    double currentUpper = upper_[iSequence];
-    double alpha = array[iPivot];
-    assert (currentSolution >= currentLower - 100.0*primalTolerance_);
-    assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
-    double thetaCoefficient;
-    double hitsLower = COIN_DBL_MAX;
-    thetaCoefficient = lowerChange[iSequence] + alpha;
-    if (thetaCoefficient > 1.0e-8)
-      hitsLower = (currentSolution - currentLower) / thetaCoefficient;
-    //if (hitsLower < 0.0) {
-    // does not hit - but should we check further
-    //   hitsLower = COIN_DBL_MAX;
-    //}
-    double hitsUpper = COIN_DBL_MAX;
-    thetaCoefficient = upperChange[iSequence] + alpha;
-    if (thetaCoefficient < -1.0e-8)
-      hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
-    //if (hitsUpper < 0.0) {
-    // does not hit - but should we check further
-    //   hitsUpper = COIN_DBL_MAX;
-    //}
-    if (CoinMin(hitsLower, hitsUpper) < theta_) {
-      theta_ = CoinMin(hitsLower, hitsUpper);
-      toLower = hitsLower < hitsUpper;
-      pivotRow_ = iPivot;
-    }
-  }
-#else
+  char * markDone = paramData.markDone;
+  const int * backwardBasic = paramData.backwardBasic;
   // first ones with alpha
   for (int i=0;i<number;i++) {
@@ -4307 +4349 @@
     // bounds will be bounds + change *theta
     double currentSolution = solution_[iSequence];
-    double currentLower = lower_[iSequence];
-    double currentUpper = upper_[iSequence];
     double alpha = array[iPivot];
-    assert (currentSolution >= currentLower - 100.0*primalTolerance_);
-    assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
-    double thetaCoefficient;
-    thetaCoefficient = lowerChange[iSequence] + alpha;
-    if (thetaCoefficient > 1.0e-8) {
+    double thetaCoefficientLower = lowerChange[iSequence] + alpha;
+    double thetaCoefficientUpper = upperChange[iSequence] + alpha;
+    if (thetaCoefficientLower > 1.0e-8) {
+      double currentLower = lower_[iSequence];
+      assert (currentSolution >= currentLower - 100.0*primalTolerance_);
       double gap=currentSolution-currentLower;
-      if (thetaCoefficient*theta_>gap) {
-        theta_ = gap/thetaCoefficient;
+      if (thetaCoefficientLower*theta_>gap) {
+        theta_ = gap/thetaCoefficientLower;
         toLower=true;
         pivotRow_=iPivot;
       }
     }
-    thetaCoefficient = upperChange[iSequence] + alpha;
-    if (thetaCoefficient < -1.0e-8) {
+    if (thetaCoefficientUpper < -1.0e-8) {
+      double currentUpper = upper_[iSequence];
+      assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
       double gap=currentSolution-currentUpper; //negative
-      if (thetaCoefficient*theta_<gap) {
-        theta_ = gap/thetaCoefficient;
+      if (thetaCoefficientUpper*theta_<gap) {
+        theta_ = gap/thetaCoefficientUpper;
         toLower=false;
         pivotRow_=iPivot;
@@ -4346 +4387 @@
           theta_ = gap/thetaCoefficient;
           toLower=true;
-          pivotRow_ = -2-iSequence;
+          pivotRow_ = backwardBasic[iSequence];
         }
       }
@@ -4364 +4405 @@
           theta_ = gap/thetaCoefficient;
           toLower=false;
-          pivotRow_ = -2-iSequence;
+          pivotRow_ = backwardBasic[iSequence];
         }
       }
     }
   }
-  if (pivotRow_<-1) {
-    // find
-    int iSequence = -pivotRow_-2;
-    for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
-      if (iSequence == pivotVariable_[iPivot]) {
-        pivotRow_=iPivot;
-        break;
-      }
-    }
-    assert (pivotRow_>=0);
-  }
-#endif
   theta_ = CoinMax(theta_,0.0);
   if (theta_>1.0e-15) {
@@ -6707 +6736 @@
   double * weights = new double [numberLook+numberColumns_];
   double * columnWeights = weights+numberLook;
+#ifndef COIN_REUSE_RANDOM
   coin_init_random_vec(columnWeights,numberColumns_);
+#else
+  for (int i=0;i<numberColumns_;i++)
+    columnWeights[i]=CoinDrand48();
+#endif
   // get row copy
   CoinPackedMatrix rowCopy = *matrix();
@@ -6961 +6995 @@
         }
       }
-      if (sum==0.0||sum>=lastSum)
+      if (sum==0.0||sum>=lastSum-1.0e-8)
         break;
       lastSum=sum;

trunk/Clp/src/ClpSimplexOther.hpp

@@ -103 +103 @@
     /// Finds best possible pivot
     double bestPivot(bool justColumns=false);
+  typedef struct {
+    double startingTheta;
+    double endingTheta;
+    double * lowerChange; // full array of lower bound changes
+    int * lowerList; // list of lower bound changes
+    double * upperChange; // full array of upper bound changes
+    int * upperList; // list of upper bound changes
+    char * markDone; // mark which ones looked at
+    int * backwardBasic; // from sequence to pivot row
+  } parametricsData;
 
 private:
@@ -113 +123 @@
          if event handler exists it may do more
      */
-     int parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
+     int parametricsLoop(parametricsData & paramData, double reportIncrement,
                          const double * changeLower, const double * changeUpper,
                          const double * changeObjective, ClpDataSave & data,
                          bool canTryQuick);
-     int parametricsLoop(double & startingTheta, double & endingTheta,
+     int parametricsLoop(parametricsData & paramData,
                          ClpDataSave & data,bool canSkipFactorization=false);
      /**  Refactorizes if necessary
@@ -137 +147 @@
          +3 max iterations
       */
-     int whileIterating(double & startingTheta, double & endingTheta, double reportIncrement,
-                        const double * changeLower, const double * changeUpper,
+     int whileIterating(parametricsData & paramData, double reportIncrement,
                         const double * changeObjective);
      /** Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none).
@@ -144 +153 @@
          type 1 bounds, 2 objective, 3 both.
      */
-     int nextTheta(int type, double maxTheta, 
-                   const double * changeLower, const double * changeUpper,
+     int nextTheta(int type, double maxTheta, parametricsData & paramData,
                    const double * changeObjective);
      /// Restores bound to original bound
@@ -197 +205 @@
          non-zero return code indicates minor problems
      */
-     int restoreFromDual(const ClpSimplex * dualProblem);
+  int restoreFromDual(const ClpSimplex * dualProblem,
+                      bool checkAccuracy=false);
      /** Does very cursory presolve.
          rhs is numberRows, whichRows is 3*numberRows and whichColumns is 2*numberColumns.