Changeset 180 for branches/pre

Timestamp:

Jul 18, 2003 3:17:33 PM (18 years ago)

Author:

forrest

Message:

new stuff

Location:

branches/pre

Files:

• ClpDualRowDantzig.cpp (modified) (2 diffs)
• ClpDualRowPivot.cpp (modified) (1 diff)
• ClpDualRowSteepest.cpp (modified) (11 diffs)
• ClpFactorization.cpp (modified) (4 diffs)
• ClpLinearObjective.cpp (modified) (2 diffs)
• ClpMatrixBase.cpp (modified) (1 diff)
• ClpMessage.cpp (modified) (2 diffs)
• ClpModel.cpp (modified) (16 diffs)
• ClpNetworkBasis.cpp (modified) (4 diffs)
• ClpNetworkMatrix.cpp (modified) (4 diffs)
• ClpNonLinearCost.cpp (modified) (5 diffs)
• ClpObjective.cpp (modified) (1 diff)
• ClpPackedMatrix.cpp (modified) (13 diffs)
• ClpPlusMinusOneMatrix.cpp (modified) (10 diffs)
• ClpPrimalColumnDantzig.cpp (modified) (1 diff)
• ClpPrimalColumnSteepest.cpp (modified) (16 diffs)
• ClpSimplex.cpp (modified) (90 diffs)
• ClpSimplexDual.cpp (modified) (93 diffs)
• ClpSimplexPrimal.cpp (modified) (38 diffs)
• ClpSimplexPrimalQuadratic.cpp (modified) (38 diffs)
• Idiot.cpp (modified) (8 diffs)
• Samples/driver.cpp (modified) (4 diffs)
• Samples/ekk.cpp (modified) (1 diff)
• Samples/ekk_interface.cpp (modified) (4 diffs)
• Samples/network.cpp (modified) (1 diff)
• Samples/useVolume.cpp (modified) (1 diff)
• Test/ClpMain.cpp (modified) (16 diffs)
• Test/MyMessageHandler.cpp (modified) (1 diff)
• Test/unitTest.cpp (modified) (19 diffs)
• include/ClpDualRowDantzig.hpp (modified) (1 diff)
• include/ClpDualRowPivot.hpp (modified) (2 diffs)
• include/ClpDualRowSteepest.hpp (modified) (4 diffs)
• include/ClpFactorization.hpp (modified) (3 diffs)
• include/ClpLinearObjective.hpp (modified) (2 diffs)
• include/ClpMatrixBase.hpp (modified) (6 diffs)
• include/ClpMessage.hpp (modified) (1 diff)
• include/ClpModel.hpp (modified) (5 diffs)
• include/ClpNetworkBasis.hpp (modified) (1 diff)
• include/ClpNetworkMatrix.hpp (modified) (3 diffs)
• include/ClpNonLinearCost.hpp (modified) (1 diff)
• include/ClpObjective.hpp (modified) (1 diff)
• include/ClpPackedMatrix.hpp (modified) (7 diffs)
• include/ClpPlusMinusOneMatrix.hpp (modified) (5 diffs)
• include/ClpPrimalColumnPivot.hpp (modified) (1 diff)
• include/ClpPrimalColumnSteepest.hpp (modified) (1 diff)
• include/ClpSimplex.hpp (modified) (16 diffs)
• include/ClpSimplexDual.hpp (modified) (2 diffs)
• include/ClpSimplexPrimal.hpp (modified) (2 diffs)
• include/ClpSimplexPrimalQuadratic.hpp (modified) (4 diffs)

branches/pre/ClpDualRowDantzig.cpp

@@ -76 +76 @@
   return chosenRow;
 }
+// Returns pivot alpha
+double
+ClpDualRowDantzig::updateWeights(CoinIndexedVector * input,
+                                  CoinIndexedVector * spare,
+                                  CoinIndexedVector * updatedColumn)
+{
+  // pivot element
+  double alpha=0.0;
+  // look at updated column
+  double * work = updatedColumn->denseVector();
+  int number = updatedColumn->getNumElements();
+  int * which = updatedColumn->getIndices();
+  int i;
+  int pivotRow = model_->pivotRow();
+
+  if (updatedColumn->packedMode()) {
+    for (i =0; i < number; i++) {
+      int iRow = which[i];
+      if (iRow==pivotRow) {
+        alpha = work[i];
+        break;
+      }
+    }
+  } else {
+    alpha = work[pivotRow];
+  }
+  return alpha;
+}
   
 /* Updates primal solution (and maybe list of candidates)
@@ -92 +120 @@
   double changeObj=0.0;
   const int * pivotVariable = model_->pivotVariable();
-  for (i=0;i<number;i++) {
-    int iRow=which[i];
-    int iPivot=pivotVariable[iRow];
-    double & value = model_->solutionAddress(iPivot);
-    double cost = model_->cost(iPivot);
-    double change = primalRatio*work[iRow];
-    value -= change;
-    changeObj -= change*cost;
-    work[iRow]=0.0;
+  if (primalUpdate->packedMode()) {
+    for (i=0;i<number;i++) {
+      int iRow=which[i];
+      int iPivot=pivotVariable[iRow];
+      double & value = model_->solutionAddress(iPivot);
+      double cost = model_->cost(iPivot);
+      double change = primalRatio*work[i];
+      value -= change;
+      changeObj -= change*cost;
+      work[i]=0.0;
+    }
+  } else {
+    for (i=0;i<number;i++) {
+      int iRow=which[i];
+      int iPivot=pivotVariable[iRow];
+      double & value = model_->solutionAddress(iPivot);
+      double cost = model_->cost(iPivot);
+      double change = primalRatio*work[iRow];
+      value -= change;
+      changeObj -= change*cost;
+      work[iRow]=0.0;
+    }
   }
   primalUpdate->setNumElements(0);

branches/pre/ClpDualRowPivot.cpp

@@ -60 +60 @@
 {
 }
-
-void 
-ClpDualRowPivot::updateWeights(CoinIndexedVector * input,
-                               CoinIndexedVector * spare,
-                               CoinIndexedVector * updatedColumn)
-{
-}
 void 
 ClpDualRowPivot::unrollWeights()

branches/pre/ClpDualRowSteepest.cpp

@@ -135 +135 @@
   // allow tolerance at least slightly bigger than standard
   tolerance = tolerance +  error;
+  // But cap
+  tolerance = min(1000.0,tolerance);
   tolerance *= tolerance; // as we are using squares
   double * solution = model_->solutionRegion();
@@ -156 +158 @@
       if (iPivot<numberColumns)
         value *= COLUMN_BIAS; // bias towards columns
+k
 #endif
       // store square in list
@@ -181 +184 @@
     number = infeasible_->getNumElements();
   }
-  for (i=0;i<number;i++) {
-    iRow = index[i];
-    double value = infeas[iRow];
-    if (value>largest*weights_[iRow]&&value>tolerance) {
-      // make last pivot row last resort choice
-      if (iRow==lastPivotRow) {
-        if (value*1.0e-10<largest*weights_[iRow]) 
-          continue;
-        else 
-          value *= 1.0e-10;
-      }
-      int iSequence = pivotVariable[iRow];
-      if (!model_->flagged(iSequence)) {
-        //#define CLP_DEBUG 1
+  if(model_->numberIterations()<model_->lastBadIteration()+200) {
+    // we can't really trust infeasibilities if there is dual error
+    double checkTolerance = 1.0e-8;
+    if (!model_->factorization()->pivots())
+      checkTolerance = 1.0e-6;
+    if (model_->largestPrimalError()>checkTolerance)
+      tolerance *= model_->largestPrimalError()/checkTolerance;
+  }
+  int numberWanted;
+  if (mode_<2)
+    numberWanted = number+1;
+  else
+    numberWanted = max(2000,number/8);
+  int iPass;
+  // Setup two passes
+  int start[4];
+  start[1]=number;
+  start[2]=0;
+  double dstart = ((double) number) * CoinDrand48();
+  start[0]=(int) dstart;
+  start[3]=start[0];
+  for (iPass=0;iPass<2;iPass++) {
+    int end = start[2*iPass+1];
+    for (i=start[2*iPass];i<end;i++) {
+      iRow = index[i];
+      double value = infeas[iRow];
+      if (value>tolerance) {
+        if (value>largest*weights_[iRow]) {
+          // make last pivot row last resort choice
+          if (iRow==lastPivotRow) {
+            if (value*1.0e-10<largest*weights_[iRow]) 
+              continue;
+            else 
+              value *= 1.0e-10;
+          }
+          int iSequence = pivotVariable[iRow];
+          if (!model_->flagged(iSequence)) {
+            //#define CLP_DEBUG 1
 #ifdef CLP_DEBUG
-        double value2=0.0;
-        if (solution[iSequence]>upper[iSequence]+tolerance) 
-          value2=solution[iSequence]-upper[iSequence];
-        else if (solution[iSequence]<lower[iSequence]-tolerance) 
-          value2=solution[iSequence]-lower[iSequence];
-        assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
-#endif
-        if (solution[iSequence]>upper[iSequence]+tolerance||
-            solution[iSequence]<lower[iSequence]-tolerance) {
-          chosenRow=iRow;
-          largest=value/weights_[iRow];
+            double value2=0.0;
+            if (solution[iSequence]>upper[iSequence]+tolerance) 
+              value2=solution[iSequence]-upper[iSequence];
+            else if (solution[iSequence]<lower[iSequence]-tolerance) 
+              value2=solution[iSequence]-lower[iSequence];
+            assert(fabs(value2*value2-infeas[iRow])<1.0e-8*min(value2*value2,infeas[iRow]));
+#endif
+            if (solution[iSequence]>upper[iSequence]+tolerance||
+                solution[iSequence]<lower[iSequence]-tolerance) {
+              chosenRow=iRow;
+              largest=value/weights_[iRow];
+            }
+          }
         }
-      }
-    }
+        numberWanted--;
+        if (!numberWanted)
+          break;
+      }
+    }
+    if (!numberWanted)
+      break;
   }
   return chosenRow;
 }
 #define TRY_NORM 1.0e-4
-// Updates weights 
-void 
+// Updates weights and returns pivot alpha
+double
 ClpDualRowSteepest::updateWeights(CoinIndexedVector * input,
                                   CoinIndexedVector * spare,
                                   CoinIndexedVector * updatedColumn)
 {
-  // clear other region
-  alternateWeights_->clear();
-  double norm = 0.0;
-  double * work = input->denseVector();
-  int number = input->getNumElements();
-  int * which = input->getIndices();
-  double * work2 = spare->denseVector();
-  int * which2 = spare->getIndices();
-  double * work3 = alternateWeights_->denseVector();
-  int * which3 = alternateWeights_->getIndices();
-  int i;
 #if CLP_DEBUG>2
   // Very expensive debug
@@ -261 +284 @@
   }
 #endif
-  for (i=0;i<number;i++) {
-    int iRow = which[i];
-    double value = work[iRow];
-    norm += value*value;
-    work2[iRow]=value;
-    which2[i]=iRow;
-  }
-  spare->setNumElements(number);
-  // ftran
-  model_->factorization()->updateColumn(alternateWeights_,spare);
-  // alternateWeights_ should still be empty
-  int pivotRow = model_->pivotRow();
+  assert (input->packedMode());
+  assert (updatedColumn->packedMode());
+  double alpha=0.0;
+  if (!model_->factorization()->networkBasis()) {
+    // clear other region
+    alternateWeights_->clear();
+    double norm = 0.0;
+    int i;
+    double * work = input->denseVector();
+    int numberNonZero = input->getNumElements();
+    int * which = input->getIndices();
+    double * work2 = spare->denseVector();
+    int * which2 = spare->getIndices();
+    // ftran
+    //permute and move indices into index array
+    //also compute norm
+    //int *regionIndex = alternateWeights_->getIndices (  );
+    const int *permute = model_->factorization()->permute();
+    //double * region = alternateWeights_->denseVector();
+    for ( i = 0; i < numberNonZero; i ++ ) {
+      int iRow = which[i];
+      double value = work[i];
+      norm += value*value;
+      iRow = permute[iRow];
+      work2[iRow] = value;
+      which2[i] = iRow;
+    }
+    spare->setNumElements ( numberNonZero );
+    model_->factorization()->updateColumn(spare,spare,true);
+    // permute back
+    numberNonZero = spare->getNumElements();
+    // alternateWeights_ should still be empty
+    int pivotRow = model_->pivotRow();
+    // could re-initialize here (could be expensive)
+    norm /= model_->alpha() * model_->alpha();
+    
+    assert(norm);
+    // pivot element
+    alpha=0.0;
+    double multiplier = 2.0 / model_->alpha();
+    // look at updated column
+    work = updatedColumn->denseVector();
+    numberNonZero = updatedColumn->getNumElements();
+    which = updatedColumn->getIndices();
+    
+    int nSave=0;
+    double * work3 = alternateWeights_->denseVector();
+    int * which3 = alternateWeights_->getIndices();
+    const int * pivotColumn = model_->factorization()->pivotColumn();
+    for (i =0; i < numberNonZero; i++) {
+      int iRow = which[i];
+      double theta = work[i];
+      if (iRow==pivotRow)
+        alpha = theta;
+      double devex = weights_[iRow];
+      work3[nSave]=devex; // save old
+      which3[nSave++]=iRow;
+      // transform to match spare
+      int jRow = pivotColumn[iRow];
+      double value = work2[jRow];
+      devex +=  theta * (theta*norm+value * multiplier);
+      if (devex < TRY_NORM) 
+        devex = TRY_NORM;
+      weights_[iRow]=devex;
+    }
+    assert (alpha);
+    alternateWeights_->setPackedMode(true);
+    alternateWeights_->setNumElements(nSave);
+    if (norm < TRY_NORM) 
+      norm = TRY_NORM;
+    weights_[pivotRow] = norm;
+    spare->clear();
 #ifdef CLP_DEBUG
-  if ( model_->logLevel (  ) >4  && 
-       fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
-    printf("on row %d, true weight %g, old %g\n",
-           pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
-#endif
-  // could re-initialize here (could be expensive)
-  norm /= model_->alpha() * model_->alpha();
-
-  assert(norm);
-  if (norm < TRY_NORM) 
-    norm = TRY_NORM;
-  if (norm != 0.) {
+    spare->checkClear();
+#endif
+  } else {
+    // clear other region
+    alternateWeights_->clear();
+    double norm = 0.0;
+    int i;
+    double * work = input->denseVector();
+    int number = input->getNumElements();
+    int * which = input->getIndices();
+    double * work2 = spare->denseVector();
+    int * which2 = spare->getIndices();
+    for (i=0;i<number;i++) {
+      int iRow = which[i];
+      double value = work[i];
+      norm += value*value;
+      work2[iRow]=value;
+      which2[i]=iRow;
+    }
+    spare->setNumElements(number);
+    // ftran
+    model_->factorization()->updateColumn(alternateWeights_,spare);
+    // alternateWeights_ should still be empty
+    int pivotRow = model_->pivotRow();
+#ifdef CLP_DEBUG
+    if ( model_->logLevel (  ) >4  && 
+         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
+      printf("on row %d, true weight %g, old %g\n",
+             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
+#endif
+    // could re-initialize here (could be expensive)
+    norm /= model_->alpha() * model_->alpha();
+    
+    assert(norm);
+    //if (norm < TRY_NORM) 
+    //norm = TRY_NORM;
+    // pivot element
+    alpha=0.0;
     double multiplier = 2.0 / model_->alpha();
     // look at updated column
@@ -291 +400 @@
     number = updatedColumn->getNumElements();
     which = updatedColumn->getIndices();
-
+    
     int nSave=0;
-
+    double * work3 = alternateWeights_->denseVector();
+    int * which3 = alternateWeights_->getIndices();
     for (i =0; i < number; i++) {
       int iRow = which[i];
-      double theta = work[iRow];
-      if (theta) {
-        double devex = weights_[iRow];
-        work3[iRow]=devex; // save old
-        which3[nSave++]=iRow;
-        double value = work2[iRow];
-        devex +=  theta * (theta*norm+value * multiplier);
-        if (devex < TRY_NORM) 
-          devex = TRY_NORM;
-        weights_[iRow]=devex;
-      }
-    }
-#ifdef CLP_DEBUG
-    assert(work3[pivotRow]&&work[pivotRow]);
-#endif
+      double theta = work[i];
+      if (iRow==pivotRow)
+        alpha = theta;
+      double devex = weights_[iRow];
+      work3[nSave]=devex; // save old
+      which3[nSave++]=iRow;
+      double value = work2[iRow];
+      devex +=  theta * (theta*norm+value * multiplier);
+      if (devex < TRY_NORM) 
+        devex = TRY_NORM;
+      weights_[iRow]=devex;
+    }
+    assert (alpha);
+    alternateWeights_->setPackedMode(true);
     alternateWeights_->setNumElements(nSave);
     if (norm < TRY_NORM) 
       norm = TRY_NORM;
     weights_[pivotRow] = norm;
-  }
-  spare->clear();
+    spare->clear();
+  }
 #ifdef CLP_DEBUG
   spare->checkClear();
 #endif
+  return alpha;
 }
   
@@ -345 +455 @@
   int numberColumns = model_->numberColumns();
 #endif
-  for (i=0;i<number;i++) {
-    int iRow=which[i];
-    int iPivot=pivotVariable[iRow];
-    double value = solution[iPivot];
-    double cost = model_->cost(iPivot);
-    double change = primalRatio*work[iRow];
-    value -= change;
-    changeObj -= change*cost;
-    solution[iPivot] = value;
-    double lower = model_->lower(iPivot);
-    double upper = model_->upper(iPivot);
-    // But if pivot row then use value of incoming
-    // Although it is safer to recompute before next selection
-    // in case something odd happens
-    if (iRow==pivotRow) {
-      iPivot = model_->sequenceIn();
-      lower = model_->lower(iPivot);
-      upper = model_->upper(iPivot);
-      value = model_->valueIncomingDual();
-    }
-    if (value<lower-tolerance) {
-      value -= lower;
-      value *= value;
-#ifdef COLUMN_BIAS 
-      if (iPivot<numberColumns)
-        value *= COLUMN_BIAS; // bias towards columns
+  if (primalUpdate->packedMode()) {
+    for (i=0;i<number;i++) {
+      int iRow=which[i];
+      int iPivot=pivotVariable[iRow];
+      double value = solution[iPivot];
+      double cost = model_->cost(iPivot);
+      double change = primalRatio*work[i];
+      work[i]=0.0;
+      value -= change;
+      changeObj -= change*cost;
+      solution[iPivot] = value;
+      double lower = model_->lower(iPivot);
+      double upper = model_->upper(iPivot);
+      // But if pivot row then use value of incoming
+      // Although it is safer to recompute before next selection
+      // in case something odd happens
+      if (iRow==pivotRow) {
+        iPivot = model_->sequenceIn();
+        lower = model_->lower(iPivot);
+        upper = model_->upper(iPivot);
+        value = model_->valueIncomingDual();
+      }
+      if (value<lower-tolerance) {
+        value -= lower;
+        value *= value;
+#ifdef COLUMN_BIAS 
+        if (iPivot<numberColumns)
+          value *= COLUMN_BIAS; // bias towards columns
 #endif
 #ifdef FIXED_BIAS 
-      if (lower==upper)
-        value *= FIXED_BIAS; // bias towards taking out fixed variables
-#endif
-      // store square in list
-      if (infeas[iRow])
-        infeas[iRow] = value; // already there
-      else
-        infeasible_->quickAdd(iRow,value);
-    } else if (value>upper+tolerance) {
-      value -= upper;
-      value *= value;
-#ifdef COLUMN_BIAS 
-      if (iPivot<numberColumns)
-        value *= COLUMN_BIAS; // bias towards columns
+        if (lower==upper)
+          value *= FIXED_BIAS; // bias towards taking out fixed variables
+#endif
+        // store square in list
+        if (infeas[iRow])
+          infeas[iRow] = value; // already there
+        else
+          infeasible_->quickAdd(iRow,value);
+      } else if (value>upper+tolerance) {
+        value -= upper;
+        value *= value;
+#ifdef COLUMN_BIAS 
+        if (iPivot<numberColumns)
+          value *= COLUMN_BIAS; // bias towards columns
 #endif
 #ifdef FIXED_BIAS 
-      if (lower==upper)
-        value *= FIXED_BIAS; // bias towards taking out fixed variables
-#endif
-      // store square in list
-      if (infeas[iRow])
-        infeas[iRow] = value; // already there
-      else
-        infeasible_->quickAdd(iRow,value);
-    } else {
-      // feasible - was it infeasible - if so set tiny
-      if (infeas[iRow])
-        infeas[iRow] = 1.0e-100;
-    }
-    work[iRow]=0.0;
+        if (lower==upper)
+          value *= FIXED_BIAS; // bias towards taking out fixed variables
+#endif
+        // store square in list
+        if (infeas[iRow])
+          infeas[iRow] = value; // already there
+        else
+          infeasible_->quickAdd(iRow,value);
+      } else {
+        // feasible - was it infeasible - if so set tiny
+        if (infeas[iRow])
+          infeas[iRow] = 1.0e-100;
+      }
+    }
+  } else {
+    for (i=0;i<number;i++) {
+      int iRow=which[i];
+      int iPivot=pivotVariable[iRow];
+      double value = solution[iPivot];
+      double cost = model_->cost(iPivot);
+      double change = primalRatio*work[iRow];
+      value -= change;
+      changeObj -= change*cost;
+      solution[iPivot] = value;
+      double lower = model_->lower(iPivot);
+      double upper = model_->upper(iPivot);
+      // But if pivot row then use value of incoming
+      // Although it is safer to recompute before next selection
+      // in case something odd happens
+      if (iRow==pivotRow) {
+        iPivot = model_->sequenceIn();
+        lower = model_->lower(iPivot);
+        upper = model_->upper(iPivot);
+        value = model_->valueIncomingDual();
+      }
+      if (value<lower-tolerance) {
+        value -= lower;
+        value *= value;
+#ifdef COLUMN_BIAS 
+        if (iPivot<numberColumns)
+          value *= COLUMN_BIAS; // bias towards columns
+#endif
+#ifdef FIXED_BIAS 
+        if (lower==upper)
+          value *= FIXED_BIAS; // bias towards taking out fixed variables
+#endif
+        // store square in list
+        if (infeas[iRow])
+          infeas[iRow] = value; // already there
+        else
+          infeasible_->quickAdd(iRow,value);
+      } else if (value>upper+tolerance) {
+        value -= upper;
+        value *= value;
+#ifdef COLUMN_BIAS 
+        if (iPivot<numberColumns)
+          value *= COLUMN_BIAS; // bias towards columns
+#endif
+#ifdef FIXED_BIAS 
+        if (lower==upper)
+          value *= FIXED_BIAS; // bias towards taking out fixed variables
+#endif
+        // store square in list
+        if (infeas[iRow])
+          infeas[iRow] = value; // already there
+        else
+          infeasible_->quickAdd(iRow,value);
+      } else {
+        // feasible - was it infeasible - if so set tiny
+        if (infeas[iRow])
+          infeas[iRow] = 1.0e-100;
+      }
+      work[iRow]=0.0;
+    }
   }
   primalUpdate->setNumElements(0);
@@ -442 +614 @@
       alternateWeights_->reserve(numberRows+
                                  model_->factorization()->maximumPivots());
-      if (!mode_||mode==5) {
+      if (!mode_||mode_==2||mode==5) {
         // initialize to 1.0 (can we do better?)
         for (i=0;i<numberRows;i++) {
@@ -455 +627 @@
         for (i=0;i<numberRows;i++) {
           double value=0.0;
-          array[i] = 1.0;
+          array[0] = 1.0;
           which[0] = i;
           alternateWeights_->setNumElements(1);
+          alternateWeights_->setPackedMode(true);
           model_->factorization()->updateColumnTranspose(temp,
                                                          alternateWeights_);
@@ -463 +636 @@
           int j;
           for (j=0;j<number;j++) {
-            int iRow=which[j];
-            value+=array[iRow]*array[iRow];
-            array[iRow]=0.0;
+            value+=array[j]*array[j];
+            array[j]=0.0;
           }
           alternateWeights_->setNumElements(0);
@@ -566 +738 @@
   int * which = alternateWeights_->getIndices();
   int i;
-  for (i=0;i<number;i++) {
-    int iRow = which[i];
-    weights_[iRow]=saved[iRow];
-    saved[iRow]=0.0;
+  if (alternateWeights_->packedMode()) {
+    for (i=0;i<number;i++) {
+      int iRow = which[i];
+      weights_[iRow]=saved[i];
+      saved[i]=0.0;
+    }
+  } else {
+    for (i=0;i<number;i++) {
+      int iRow = which[i];
+      weights_[iRow]=saved[iRow];
+      saved[iRow]=0.0;
+    }
   }
   alternateWeights_->setNumElements(0);
@@ -598 +778 @@
   state_ =-1;
 }
+// Returns true if would not find any row
+bool 
+ClpDualRowSteepest::looksOptimal() const
+{
+  int iRow;
+  const int * pivotVariable = model_->pivotVariable();
+  double tolerance=model_->currentPrimalTolerance();
+  // we can't really trust infeasibilities if there is primal error
+  // this coding has to mimic coding in checkPrimalSolution
+  double error = min(1.0e-3,model_->largestPrimalError());
+  // allow tolerance at least slightly bigger than standard
+  tolerance = tolerance +  error;
+  // But cap
+  tolerance = min(1000.0,tolerance);
+  int numberRows = model_->numberRows();
+  int numberInfeasible=0;
+  for (iRow=0;iRow<numberRows;iRow++) {
+    int iPivot=pivotVariable[iRow];
+    double value = model_->solution(iPivot);
+    double lower = model_->lower(iPivot);
+    double upper = model_->upper(iPivot);
+    if (value<lower-tolerance) {
+      numberInfeasible++;
+    } else if (value>upper+tolerance) {
+      numberInfeasible++;
+    }
+  }
+  return (numberInfeasible==0);
+}
 

branches/pre/ClpFactorization.cpp

@@ -73 +73 @@
 }
 int 
-ClpFactorization::factorize ( const ClpSimplex * model,
-                              const ClpMatrixBase * matrix, 
-                              int numberRows, int numberColumns,
-                              int rowIsBasic[], int columnIsBasic[] , 
-                              double areaFactor )
-{
+ClpFactorization::factorize ( ClpSimplex * model,
+                              int solveType, bool valuesPass)
+{
+  const ClpMatrixBase * matrix = model->clpMatrix(); 
+  int numberRows = model->numberRows();
+  int numberColumns = model->numberColumns();
+  // If too many compressions increase area
+  if (numberPivots_>1&&numberCompressions_*10 > numberPivots_+10) {
+    areaFactor_ *= 1.1;
+  }
   if (!networkBasis_||doCheck) {
-    // see if matrix a network
-    ClpNetworkMatrix* networkMatrix =
-      dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
-    // If network - still allow ordinary factorization first time for laziness
-    int saveMaximumPivots = maximumPivots();
-    delete networkBasis_;
-    networkBasis_ = NULL;
-    if (networkMatrix&&!doCheck)
-      maximumPivots(1);
-    // maybe for speed will be better to leave as many regions as possible
-    gutsOfDestructor();
-    gutsOfInitialize(2);
-    if (areaFactor)
-      areaFactor_ = areaFactor;
-    int numberBasic = 0;
-    CoinBigIndex numberElements=0;
-    int numberRowBasic=0;
-    
-    // compute how much in basis
-    
-    int i;
-    
-    for (i=0;i<numberRows;i++) {
-      if (rowIsBasic[i]>=0)
-        numberRowBasic++;
-    }
-    
-    numberBasic = numberRowBasic;
-    for (i=0;i<numberColumns;i++) {
-      if (columnIsBasic[i]>=0) {
-        numberBasic++;
+    status_=-99;
+    int * pivotVariable = model->pivotVariable();
+    //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
+    while (status_<-98) {
+      
+      int i;
+      int numberBasic=0;
+      int numberRowBasic;
+      for (i=0;i<numberRows;i++) {
+        if (model->getRowStatus(i) == ClpSimplex::basic) 
+          pivotVariable[numberBasic++]=i;
       }
-    }
-    numberElements += matrix->numberInBasis(columnIsBasic);
-    if ( numberBasic > numberRows ) {
-    return -2; // say too many in basis
-    }
-    numberElements = 3 * numberBasic + 3 * numberElements + 10000;
-    getAreas ( numberRows, numberBasic, numberElements,
-               2 * numberElements );
-    //fill
-    //copy
-    numberBasic=0;
-    numberElements=0;
-    for (i=0;i<numberRows;i++) {
-      if (rowIsBasic[i]>=0) {
-        indexRowU_[numberElements]=i;
-        indexColumnU_[numberElements]=numberBasic;
-        elementU_[numberElements++]=slackValue_;
-        numberBasic++;
+      numberRowBasic=numberBasic;
+      for (i=0;i<numberColumns;i++) {
+        if (model->getColumnStatus(i) == ClpSimplex::basic) 
+          pivotVariable[numberBasic++]=i;
       }
-    }
-    numberElements +=matrix->fillBasis(model, columnIsBasic, numberBasic, 
-                                       indexRowU_+numberElements, 
-                                       indexColumnU_+numberElements,
-                                       elementU_+numberElements);
-    lengthU_ = numberElements;
-    
-    preProcess ( 0 );
-    factor (  );
-    numberBasic=0;
-    if (status_ == 0) {
-      int * permuteBack = permuteBack_;
-      int * back = pivotColumnBack_;
-      for (i=0;i<numberRows;i++) {
-        if (rowIsBasic[i]>=0) {
-          rowIsBasic[i]=permuteBack[back[numberBasic++]];
+      assert (numberBasic<=numberRows);
+      // see if matrix a network
+      ClpNetworkMatrix* networkMatrix =
+        dynamic_cast< ClpNetworkMatrix*>(model->clpMatrix());
+      // If network - still allow ordinary factorization first time for laziness
+      int saveMaximumPivots = maximumPivots();
+      delete networkBasis_;
+      networkBasis_ = NULL;
+      if (networkMatrix&&!doCheck)
+        maximumPivots(1);
+      while (status_==-99) {
+        // maybe for speed will be better to leave as many regions as possible
+        gutsOfDestructor();
+        gutsOfInitialize(2);
+        CoinBigIndex numberElements=numberRowBasic;
+
+        // compute how much in basis
+
+        int i;
+
+        numberElements +=matrix->fillBasis(model,
+                                           pivotVariable+numberRowBasic, 
+                                           numberRowBasic,
+                                           numberBasic-numberRowBasic,
+                                           NULL,NULL,NULL);
+
+        numberElements = 3 * numberBasic + 3 * numberElements + 10000;
+        getAreas ( numberRows, numberBasic, numberElements,
+                   2 * numberElements );
+        //fill
+        //copy
+        numberElements=numberRowBasic;
+        for (i=0;i<numberRowBasic;i++) {
+          int iRow = pivotVariable[i];
+          indexRowU_[i]=iRow;
+          indexColumnU_[i]=i;
+          elementU_[i]=slackValue_;
+        }
+        numberElements +=matrix->fillBasis(model, 
+                                           pivotVariable+numberRowBasic, 
+                                           numberRowBasic, 
+                                           numberBasic-numberRowBasic,
+                                           indexRowU_+numberElements, 
+                                           indexColumnU_+numberElements,
+                                           elementU_+numberElements);
+        lengthU_ = numberElements;
+
+        preProcess ( 0 );
+        factor (  );
+        if (status_==-99) {
+          // get more memory
+          areaFactor(2.0*areaFactor());
         }
       }
-      for (i=0;i<numberColumns;i++) {
-        if (columnIsBasic[i]>=0) {
-          columnIsBasic[i]=permuteBack[back[numberBasic++]];
-        }
-      }
-      if (increasingRows_>1) {
+      // If we get here status is 0 or -1
+      if (status_ == 0) {
+        int * permuteBack = permuteBack_;
+        int * back = pivotColumnBack_;
+        int * pivotTemp = pivotColumn_;
+        ClpDisjointCopyN ( pivotVariable, numberRows_ , pivotTemp  );
+        // Redo pivot order
+        for (i=0;i<numberRowBasic;i++) {
+          int k = pivotTemp[i];
+          // so rowIsBasic[k] would be permuteBack[back[i]]
+          pivotVariable[permuteBack[back[i]]]=k+numberColumns;
+        }
+        for (;i<numberRows;i++) {
+          int k = pivotTemp[i];
+          // so rowIsBasic[k] would be permuteBack[back[i]]
+          pivotVariable[permuteBack[back[i]]]=k;
+        }
         // Set up permutation vector
-        if (increasingRows_<3) {
-          // these arrays start off as copies of permute
-          // (and we could use permute_ instead of pivotColumn (not back though))
-          ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
-          ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
-        }
-      } else {
-        // Set up permutation vector
-        // (we could use permute_ instead of pivotColumn (not back though))
-        for (i=0;i<numberRows_;i++) {
-          int k=pivotColumn_[i];
-          pivotColumn_[i]=pivotColumnBack_[i];
-          pivotColumnBack_[i]=k;
-        }
-      }
-    } else if (status_ == -1) {
-      // mark as basic or non basic
-      for (i=0;i<numberRows_;i++) {
-        if (rowIsBasic[i]>=0) {
-          if (pivotColumn_[numberBasic]>=0) 
-            rowIsBasic[i]=pivotColumn_[numberBasic];
-          else
-            rowIsBasic[i]=-1;
-          numberBasic++;
-        }
-      }
-      for (i=0;i<numberColumns;i++) {
-        if (columnIsBasic[i]>=0) {
-          if (pivotColumn_[numberBasic]>=0) 
-            columnIsBasic[i]=pivotColumn_[numberBasic];
-          else
-            columnIsBasic[i]=-1;
-          numberBasic++;
-        }
-      }
-    }
-    if (networkMatrix) {
-      maximumPivots(saveMaximumPivots);
-      if (!status_) {
-        // create network factorization
-        if (doCheck)
-          delete networkBasis_; // temp
-        networkBasis_ = new ClpNetworkBasis(model,numberRows_,
-                                            pivotRegion_,
-                                            permuteBack_,
-                                            startColumnU_,
-                                            numberInColumn_,
-                                            indexRowU_,
-                                            elementU_);
-        // kill off arrays in ordinary factorization
-        if (!doCheck)
-          gutsOfDestructor();
-      }
-    }
-    if (!status_&&!networkBasis_) {
-      // See if worth going sparse and when
-      if (numberFtranCounts_>100) {
-        ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
-        ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
-        ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
-        if (btranCountInput_) {
-          btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
-          btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
-          btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
+        // these arrays start off as copies of permute
+        // (and we could use permute_ instead of pivotColumn (not back though))
+        ClpDisjointCopyN ( permute_, numberRows_ , pivotColumn_  );
+        ClpDisjointCopyN ( permuteBack_, numberRows_ , pivotColumnBack_  );
+        if (networkMatrix) {
+          maximumPivots(saveMaximumPivots);
+          // create network factorization
+          if (doCheck)
+            delete networkBasis_; // temp
+          networkBasis_ = new ClpNetworkBasis(model,numberRows_,
+                                              pivotRegion_,
+                                              permuteBack_,
+                                              startColumnU_,
+                                              numberInColumn_,
+                                              indexRowU_,
+                                              elementU_);
+          // kill off arrays in ordinary factorization
+          if (!doCheck) {
+            gutsOfDestructor();
+#if 0
+            // but put back permute arrays so odd things will work
+            int numberRows = model->numberRows();
+            pivotColumnBack_ = new int [numberRows];
+            permute_ = new int [numberRows];
+            int i;
+            for (i=0;i<numberRows;i++) {
+              pivotColumnBack_[i]=i;
+              permute_[i]=i;
+            }
+#endif
+          }
         } else {
-          // odd - we have not done any btrans
-          btranAverageAfterU_ = 1.0;
-          btranAverageAfterR_ = 1.0;
-          btranAverageAfterL_ = 1.0;
-        }
-      }
-      // scale back
-
-      ftranCountInput_ *= 0.8;
-      ftranCountAfterL_ *= 0.8;
-      ftranCountAfterR_ *= 0.8;
-      ftranCountAfterU_ *= 0.8;
-      btranCountInput_ *= 0.8;
-      btranCountAfterU_ *= 0.8;
-      btranCountAfterR_ *= 0.8;
-      btranCountAfterL_ *= 0.8;
+          // See if worth going sparse and when
+          if (numberFtranCounts_>100) {
+            ftranAverageAfterL_ = max(ftranCountAfterL_/ftranCountInput_,1.0);
+            ftranAverageAfterR_ = max(ftranCountAfterR_/ftranCountAfterL_,1.0);
+            ftranAverageAfterU_ = max(ftranCountAfterU_/ftranCountAfterR_,1.0);
+            if (btranCountInput_) {
+              btranAverageAfterU_ = max(btranCountAfterU_/btranCountInput_,1.0);
+              btranAverageAfterR_ = max(btranCountAfterR_/btranCountAfterU_,1.0);
+              btranAverageAfterL_ = max(btranCountAfterL_/btranCountAfterR_,1.0);
+            } else {
+              // odd - we have not done any btrans
+              btranAverageAfterU_ = 1.0;
+              btranAverageAfterR_ = 1.0;
+              btranAverageAfterL_ = 1.0;
+            }
+          }
+          // scale back
+          
+          ftranCountInput_ *= 0.8;
+          ftranCountAfterL_ *= 0.8;
+          ftranCountAfterR_ *= 0.8;
+          ftranCountAfterU_ *= 0.8;
+          btranCountInput_ *= 0.8;
+          btranCountAfterU_ *= 0.8;
+          btranCountAfterR_ *= 0.8;
+          btranCountAfterL_ *= 0.8;
+        }
+      } else if (status_ == -1&&(solveType==0||solveType==2)) {
+        // This needs redoing as it was merged coding - does not need array
+        int numberTotal = numberRows+numberColumns;
+        int * isBasic = new int [numberTotal];
+        int * rowIsBasic = isBasic+numberColumns;
+        int * columnIsBasic = isBasic;
+        for (i=0;i<numberTotal;i++) 
+          isBasic[i]=-1;
+        for (i=0;i<numberRowBasic;i++) {
+          int iRow = pivotVariable[i];
+          rowIsBasic[iRow]=1;
+        }
+        for (;i<numberBasic;i++) {
+          int iColumn = pivotVariable[i];
+          columnIsBasic[iColumn]=1;
+        }
+        numberBasic=0;
+        for (i=0;i<numberRows;i++) 
+          pivotVariable[i]=-1;
+        // mark as basic or non basic
+        for (i=0;i<numberRows;i++) {
+          if (rowIsBasic[i]>=0) {
+            if (pivotColumn_[numberBasic]>=0) 
+              rowIsBasic[i]=pivotColumn_[numberBasic];
+            else
+              rowIsBasic[i]=-1;
+            numberBasic++;
+          }
+        }
+        for (i=0;i<numberColumns;i++) {
+          if (columnIsBasic[i]>=0) {
+            if (pivotColumn_[numberBasic]>=0) 
+              columnIsBasic[i]=pivotColumn_[numberBasic];
+            else
+              columnIsBasic[i]=-1;
+            numberBasic++;
+          }
+        }
+        // leave pivotVariable in useful form for cleaning basis
+        int * pivotVariable = model->pivotVariable();
+        for (i=0;i<numberRows;i++) {
+          pivotVariable[i]=-1;
+        }
+        for (i=0;i<numberRows;i++) {
+          if (model->getRowStatus(i) == ClpSimplex::basic) {
+            int iPivot = rowIsBasic[i];
+            if (iPivot>=0) 
+              pivotVariable[iPivot]=i+numberColumns;
+          }
+        }
+        for (i=0;i<numberColumns;i++) {
+          if (model->getColumnStatus(i) == ClpSimplex::basic) {
+            int iPivot = columnIsBasic[i];
+            if (iPivot>=0) 
+              pivotVariable[iPivot]=i;
+          }
+        }
+        delete [] isBasic;
+        double * columnLower = model->lowerRegion();
+        double * columnUpper = model->upperRegion();
+        double * columnActivity = model->solutionRegion();
+        double * rowLower = model->lowerRegion(0);
+        double * rowUpper = model->upperRegion(0);
+        double * rowActivity = model->solutionRegion(0);
+        //redo basis - first take ALL columns out
+        int iColumn;
+        double largeValue = model->largeValue();
+        for (iColumn=0;iColumn<numberColumns;iColumn++) {
+          if (model->getColumnStatus(iColumn)==ClpSimplex::basic) {
+            // take out
+            if (!valuesPass) {
+              double lower=columnLower[iColumn];
+              double upper=columnUpper[iColumn];
+              double value=columnActivity[iColumn];
+              if (lower>-largeValue||upper<largeValue) {
+                if (fabs(value-lower)<fabs(value-upper)) {
+                  model->setColumnStatus(iColumn,ClpSimplex::atLowerBound);
+                  columnActivity[iColumn]=lower;
+                } else {
+                  model->setColumnStatus(iColumn,ClpSimplex::atUpperBound);
+                  columnActivity[iColumn]=upper;
+                }
+              } else {
+                model->setColumnStatus(iColumn,ClpSimplex::isFree);
+              }
+            } else {
+              model->setColumnStatus(iColumn,ClpSimplex::superBasic);
+            }
+          }
+        }
+        int iRow;
+        for (iRow=0;iRow<numberRows;iRow++) {
+          int iSequence=pivotVariable[iRow];
+          if (iSequence>=0) {
+            // basic
+            if (iSequence>=numberColumns) {
+              // slack in - leave
+              //assert (iSequence-numberColumns==iRow);
+            } else {
+              // put back structural
+              model->setColumnStatus(iSequence,ClpSimplex::basic);
+            }
+          } else {
+            // put in slack
+            model->setRowStatus(iRow,ClpSimplex::basic);
+          }
+        }
+        // signal repeat
+        status_=-99;
+        // set fixed if they are
+        for (iRow=0;iRow<numberRows;iRow++) {
+          if (model->getRowStatus(iRow)!=ClpSimplex::basic ) {
+            if (rowLower[iRow]==rowUpper[iRow]) {
+              rowActivity[iRow]=rowLower[iRow];
+              model->setRowStatus(iRow,ClpSimplex::isFixed);
+            }
+          }
+        }
+        for (iColumn=0;iColumn<numberColumns;iColumn++) {
+          if (model->getColumnStatus(iColumn)!=ClpSimplex::basic ) {
+            if (columnLower[iColumn]==columnUpper[iColumn]) {
+              columnActivity[iColumn]=columnLower[iColumn];
+              model->setColumnStatus(iColumn,ClpSimplex::isFixed);
+            }
+          }
+        }
+      } 
     }
   } else {
@@ -285 +407 @@
    region1 starts as zero and is zero at end */
 int 
-ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
-                                 CoinIndexedVector * regionSparse2,
-                                 bool FTUpdate) 
+ClpFactorization::updateColumnFT ( CoinIndexedVector * regionSparse,
+                                   CoinIndexedVector * regionSparse2)
 {
 #ifdef CLP_DEBUG
@@ -294 +415 @@
   if (!networkBasis_) {
     collectStatistics_ = true;
-    return CoinFactorization::updateColumn(regionSparse,
-                                           regionSparse2,
-                                           FTUpdate);
+    int returnValue= CoinFactorization::updateColumnFT(regionSparse,
+                                                       regionSparse2);
     collectStatistics_ = false;
+    return returnValue;
   } else {
 #ifdef CHECK_NETWORK
     CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
-    int returnCode = CoinFactorization::updateColumn(regionSparse,
-                                                     regionSparse2, FTUpdate);
+    int returnCode = CoinFactorization::updateColumnFT(regionSparse,
+                                                       regionSparse2);
     networkBasis_->updateColumn(regionSparse,save);
     int i;
@@ -315 +436 @@
     return returnCode;
 #else
-    return networkBasis_->updateColumn(regionSparse,regionSparse2);
-#endif
-  }
-}
-/* Updates one column (FTRAN) to/from array 
+    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
+    return 1;
+#endif
+  }
+}
+/* Updates one column (FTRAN) from region2
    number returned is negative if no room
-   ** For large problems you should ALWAYS know where the nonzeros
-   are, so please try and migrate to previous method after you
-   have got code working using this simple method - thank you!
-   (the only exception is if you know input is dense e.g. rhs)
-   region starts as zero and is zero at end */
+   region1 starts as zero and is zero at end */
 int 
 ClpFactorization::updateColumn ( CoinIndexedVector * regionSparse,
-                        double array[] ) const
-{
+                                 CoinIndexedVector * regionSparse2,
+                                 bool noPermute) const
+{
+#ifdef CLP_DEBUG
+  regionSparse->checkClear();
+#endif
   if (!networkBasis_) {
-    return CoinFactorization::updateColumn(regionSparse,
-                                           array);
+    collectStatistics_ = true;
+    int returnValue= CoinFactorization::updateColumn(regionSparse,
+                                                     regionSparse2,
+                                                     noPermute);
+    collectStatistics_ = false;
+    return returnValue;
   } else {
 #ifdef CHECK_NETWORK
-    double * save = new double [numberRows_+1];
-    memcpy(save,array,(numberRows_+1)*sizeof(double));
+    CoinIndexedVector * save = new CoinIndexedVector(*regionSparse2);
     int returnCode = CoinFactorization::updateColumn(regionSparse,
-                                                     array);
-    networkBasis_->updateColumn(regionSparse, save);
+                                                     regionSparse2,
+                                                     noPermute);
+    networkBasis_->updateColumn(regionSparse,save);
     int i;
-    for (i=0;i<numberRows_;i++)
-      assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
-    delete [] save;
+    double * array = regionSparse2->denseVector();
+    double * array2 = save->denseVector();
+    for (i=0;i<numberRows_;i++) {
+      double value1 = array[i];
+      double value2 = array2[i];
+      assert (value1==value2);
+    }
+    delete save;
     return returnCode;
 #else
-    return networkBasis_->updateColumn(regionSparse, array);
-#endif
-  }
-}
-/* Updates one column transpose (BTRAN)
-   For large problems you should ALWAYS know where the nonzeros
-   are, so please try and migrate to previous method after you
-   have got code working using this simple method - thank you!
-   (the only exception is if you know input is dense e.g. dense objective)
-   returns number of nonzeros */
-int 
-ClpFactorization::updateColumnTranspose ( CoinIndexedVector * regionSparse,
-                                          double array[] ) const
-{
-  if (!networkBasis_) {
-    return CoinFactorization::updateColumnTranspose(regionSparse,
-                                                    array);
-  } else {
-#ifdef CHECK_NETWORK
-    double * save = new double [numberRows_+1];
-    memcpy(save,array,(numberRows_+1)*sizeof(double));
-    int returnCode = CoinFactorization::updateColumnTranspose(regionSparse,
-                                                              array);
-    networkBasis_->updateColumnTranspose(regionSparse, save);
-    int i;
-    for (i=0;i<numberRows_;i++)
-      assert (fabs(save[i]-array[i])<1.0e-8*(1.0+fabs(array[i])));
-    delete [] save;
-    return returnCode;
-#else
-    return networkBasis_->updateColumnTranspose(regionSparse, array);
+    networkBasis_->updateColumn(regionSparse,regionSparse2,-1);
+    return 1;
 #endif
   }

branches/pre/ClpLinearObjective.cpp

@@ -53 +53 @@
   }
 }
+/* Subset constructor.  Duplicates are allowed
+   and order is as given.
+*/
+ClpLinearObjective::ClpLinearObjective (const ClpLinearObjective &rhs,
+                                        int numberColumns, 
+                                        const int * whichColumn) 
+: ClpObjective(rhs)
+{
+  objective_=NULL;
+  numberColumns_=0;
+  if (numberColumns>0) {
+    // check valid lists
+    int numberBad=0;
+    int i;
+    for (i=0;i<numberColumns;i++)
+      if (whichColumn[i]<0||whichColumn[i]>=rhs.numberColumns_)
+        numberBad++;
+    if (numberBad)
+      throw CoinError("bad column list", "subset constructor", 
+                      "ClpLinearObjective");
+    numberColumns_ = numberColumns;
+    objective_ = new double[numberColumns_];
+    for (i=0;i<numberColumns_;i++) 
+      objective_[i]=rhs.objective_[whichColumn[i]];
+  }
+}
+  
 
 //-------------------------------------------------------------------
@@ -95 +122 @@
 {
   return new ClpLinearObjective(*this);
+}
+/* Subset clone.  Duplicates are allowed
+   and order is as given.
+*/
+ClpObjective * 
+ClpLinearObjective::subsetClone (int numberColumns, 
+                           const int * whichColumns) const
+{
+  return new ClpLinearObjective(*this, numberColumns, whichColumns);
 }
 // Resize objective

branches/pre/ClpMatrixBase.cpp

@@ -69 +69 @@
   abort();
 }
+/* Subset clone (without gaps).  Duplicates are allowed
+   and order is as given.
+   Derived classes need not provide this as it may not always make
+   sense */
+ClpMatrixBase * 
+ClpMatrixBase::subsetClone (
+                            int numberRows, const int * whichRows,
+                            int numberColumns, const int * whichColumns) const
+ 
 
+{
+  std::cerr<<"subsetClone not supported - ClpMatrixBase"<<std::endl;
+  abort();
+}

branches/pre/ClpMessage.cpp

@@ -26 +26 @@
   {CLP_SIMPLEX_FLAG,10,3,"Flagging variable %c%d"},
   {CLP_SIMPLEX_GIVINGUP,11,2,"Stopping as close enough"},
-  {CLP_DUAL_CHECKB,12,3,"Looking infeasible checking bounds with %g"},
+  {CLP_DUAL_CHECKB,12,2,"New dual bound of %g"},
   {CLP_DUAL_ORIGINAL,13,3,"Going back to original objective"},
   {CLP_SIMPLEX_PERTURB,14,1,"Perturbing problem by %g"},
@@ -46 +46 @@
   {CLP_DUAL_CHECK,106,4,"Btran alpha %g, ftran alpha %g"},
   {CLP_PRIMAL_DJ,107,4,"Btran dj %g, ftran dj %g"},
-  {CLP_PRESOLVE_COLINFEAS,501,2,"Problem is infeasible due to column %d, %g %g"},
-  {CLP_PRESOLVE_ROWINFEAS,502,2,"Problem is infeasible due to row %d, %g %g"},
-  {CLP_PRESOLVE_COLUMNBOUNDA,503,2,"Problem looks unbounded above due to column %d, %g %g"},
-  {CLP_PRESOLVE_COLUMNBOUNDB,504,2,"Problem looks unbounded below due to column %d, %g %g"},
-  {CLP_PRESOLVE_NONOPTIMAL,505,1,"Problem not optimal, resolve after postsolve"},
-  {CLP_PRESOLVE_STATS,506,1,"Presolve %d (%d) rows, %d (%d) columns and %d (%d) elements"},
-  {CLP_PRESOLVE_INFEAS,507,0,"Presolve determined that the problem was infeasible with tolerance of %g"},
-  {CLP_PRESOLVE_UNBOUND,508,0,"Presolve thinks problem is unbounded"},
-  {CLP_PRESOLVE_INFEASUNBOUND,509,0,"Presolve thinks problem is infeasible AND unbounded???"},
-  {CLP_PRESOLVE_INTEGERMODS,510,1,"Presolve is modifying %d integer bounds and re-presolving"},
-  {CLP_PRESOLVE_POSTSOLVE,511,0,"After Postsolve, objective %g, infeasibilities - dual %g (%d), primal %g (%d)"},
-  {CLP_PRESOLVE_NEEDS_CLEANING,512,0,"Presolved model was optimal, full model needs cleaning up"},
   {CLP_PACKEDSCALE_INITIAL,1001,2,"Initial range of elements is %g to %g"},
   {CLP_PACKEDSCALE_WHILE,1002,3,"Range of elements is %g to %g"},

branches/pre/ClpModel.cpp

@@ -9 +9 @@
 #include <iostream>
 
-#include <time.h>
 
 #include "CoinPragma.hpp"
-#ifndef _MSC_VER
-#include <sys/times.h>
-#include <sys/resource.h>
-#include <unistd.h>
-#endif
 
 #include "CoinHelperFunctions.hpp"
@@ -27 +21 @@
 #include "ClpLinearObjective.hpp"
 
-static double cpuTime()
-{
-  double cpu_temp;
-#if defined(_MSC_VER)
-  unsigned int ticksnow;        /* clock_t is same as int */
-  
-  ticksnow = (unsigned int)clock();
-  
-  cpu_temp = (double)((double)ticksnow/CLOCKS_PER_SEC);
-#else
-  struct rusage usage;
-  getrusage(RUSAGE_SELF,&usage);
-  cpu_temp = usage.ru_utime.tv_sec;
-  cpu_temp += 1.0e-6*((double) usage.ru_utime.tv_usec);
-#endif
-  return cpu_temp;
-}
 //#############################################################################
 
@@ -339 +316 @@
   strParam_[ClpProbName] = rhs.strParam_[ClpProbName];
 
+  optimizationDirection_ = rhs.optimizationDirection_;
   objectiveValue_=rhs.objectiveValue_;
   smallElement_ = rhs.smallElement_;
@@ -344 +322 @@
   solveType_ = rhs.solveType_;
   problemStatus_ = rhs.problemStatus_;
+  numberRows_ = rhs.numberRows_;
+  numberColumns_ = rhs.numberColumns_;
   if (trueCopy) {
     lengthNames_ = rhs.lengthNames_;
@@ -527 +507 @@
   case ClpMaxSeconds: 
     if(value>=0)
-      value += cpuTime();
+      value += CoinCpuTime();
     else
       value = -1.0;
@@ -733 +713 @@
                               number, which, newSize);
   matrix_->deleteRows(number,which);
+  //matrix_->removeGaps();
   // status
   if (status_) {
@@ -770 +751 @@
                               number, which, newSize);
   matrix_->deleteCols(number,which);
+  //matrix_->removeGaps();
   if (quadraticObjective_) {
     quadraticObjective_->deleteCols(number,which);
@@ -787 +769 @@
     status_ = temp;
   }
+  integerType_ = deleteChar(integerType_,numberColumns_,
+                            number, which, newSize,true);
   numberColumns_=newSize;
   // set state back to unknown
@@ -796 +780 @@
   rowNames_ = std::vector<std::string> ();
   columnNames_ = std::vector<std::string> ();
-  integerType_ = deleteChar(integerType_,numberColumns_,
-                            number, which, newSize,true);
 }
 // Add rows
@@ -815 +797 @@
       rows[iRow] = 
         new CoinPackedVector(rowStarts[iRow+1]-iStart,
+                             columns+iStart,elements+iStart);
+    }
+    addRows(number, rowLower, rowUpper,
+            rows);
+    for (iRow=0;iRow<number;iRow++) 
+      delete rows[iRow];
+    delete [] rows;
+  }
+}
+// Add rows
+void 
+ClpModel::addRows(int number, const double * rowLower, 
+                  const double * rowUpper,
+                  const int * rowStarts, 
+                  const int * rowLengths, const int * columns,
+                  const double * elements)
+{
+  // Create a list of CoinPackedVectors
+  if (number) {
+    CoinPackedVectorBase ** rows=
+      new CoinPackedVectorBase * [number];
+    int iRow;
+    for (iRow=0;iRow<number;iRow++) {
+      int iStart = rowStarts[iRow];
+      rows[iRow] = 
+        new CoinPackedVector(rowLengths[iRow],
                              columns+iStart,elements+iStart);
     }
@@ -886 +894 @@
       columns[iColumn] = 
         new CoinPackedVector(columnStarts[iColumn+1]-iStart,
+                             rows+iStart,elements+iStart);
+    }
+    addColumns(number, columnLower, columnUpper,
+               objIn, columns);
+    for (iColumn=0;iColumn<number;iColumn++) 
+      delete columns[iColumn];
+    delete [] columns;
+
+  }
+}
+// Add columns
+void 
+ClpModel::addColumns(int number, const double * columnLower, 
+                     const double * columnUpper,
+                     const double * objIn,
+                     const int * columnStarts, 
+                     const int * columnLengths, const int * rows,
+                     const double * elements)
+{
+  // Create a list of CoinPackedVectors
+  if (number) {
+    CoinPackedVectorBase ** columns=
+      new CoinPackedVectorBase * [number];
+    int iColumn;
+    for (iColumn=0;iColumn<number;iColumn++) {
+      int iStart = columnStarts[iColumn];
+      columns[iColumn] = 
+        new CoinPackedVector(columnLengths[iColumn],
                              rows+iStart,elements+iStart);
     }
@@ -979 +1015 @@
 {
   if(value>=0)
-    dblParam_[ClpMaxSeconds]=value+cpuTime();
+    dblParam_[ClpMaxSeconds]=value+CoinCpuTime();
   else
     dblParam_[ClpMaxSeconds]=-1.0;
@@ -989 +1025 @@
   bool hitMax= (numberIterations_>=maximumIterations());
   if (dblParam_[ClpMaxSeconds]>=0.0&&!hitMax)
-    hitMax = (cpuTime()>=dblParam_[ClpMaxSeconds]);
+    hitMax = (CoinCpuTime()>=dblParam_[ClpMaxSeconds]);
   return hitMax;
 }
@@ -1030 +1066 @@
   }
   CoinMpsIO m;
-  double time1 = cpuTime(),time2;
+  double time1 = CoinCpuTime(),time2;
   int status=m.readMps(fileName,"");
   if (!status||ignoreErrors) {
@@ -1068 +1104 @@
     }
     setDblParam(ClpObjOffset,m.objectiveOffset());
-    time2 = cpuTime();
+    time2 = CoinCpuTime();
     handler_->message(CLP_IMPORT_RESULT,messages_)
       <<fileName
@@ -1188 +1224 @@
   delete quadraticObjective_;
   quadraticObjective_ = NULL;
+}
+// Returns resized array and updates size
+double * whichDouble(double * array , int number, const int * which)
+{
+  double * newArray=NULL;
+  if (array&&number) {
+    int i ;
+    newArray = new double[number];
+    for (i=0;i<number;i++) 
+      newArray[i]=array[which[i]];
+  }
+  return newArray;
+}
+char * whichChar(char * array , int number, const int * which)
+{
+  char * newArray=NULL;
+  if (array&&number) {
+    int i ;
+    newArray = new char[number];
+    for (i=0;i<number;i++) 
+      newArray[i]=array[which[i]];
+  }
+  return newArray;
+}
+unsigned char * whichUnsignedChar(unsigned char * array , 
+                                  int number, const int * which)
+{
+  unsigned char * newArray=NULL;
+  if (array&&number) {
+    int i ;
+    newArray = new unsigned char[number];
+    for (i=0;i<number;i++) 
+      newArray[i]=array[which[i]];
+  }
+  return newArray;
+}
+// Subproblem constructor
+ClpModel::ClpModel ( const ClpModel * rhs,
+                     int numberRows, const int * whichRow,
+                     int numberColumns, const int * whichColumn,
+                     bool dropNames, bool dropIntegers)
+{
+#if 0
+  // Could be recoded to be faster and take less memory
+  // and to allow re-ordering and duplicates etc
+  gutsOfCopy(*rhs,true);
+  int numberRowsWhole = rhs->numberRows();
+  int * delRow = new int[numberRowsWhole];
+  memset(delRow,0,numberRowsWhole*sizeof(int));
+  int i;
+  for (i=0;i<numberRows;i++) {
+    int iRow=whichRow[i];
+    assert (iRow>=0&&iRow<numberRowsWhole);
+    delRow[iRow]=1;
+  }
+  numberRows=0;
+  for (i=0;i<numberRowsWhole;i++) {
+    if (delRow[i]==0)
+      delRow[numberRows++]=i;
+  }
+  deleteRows(numberRows,delRow);
+  delete [] delRow;
+  int numberColumnsWhole = rhs->numberColumns();
+  int * delColumn = new int[numberColumnsWhole];
+  memset(delColumn,0,numberColumnsWhole*sizeof(int));
+  for (i=0;i<numberColumns;i++) {
+    int iColumn=whichColumn[i];
+    assert (iColumn>=0&&iColumn<numberColumnsWhole);
+    delColumn[iColumn]=1;
+  }
+  numberColumns=0;
+  for (i=0;i<numberColumnsWhole;i++) {
+    if (delColumn[i]==0)
+      delColumn[numberColumns++]=i;
+  }
+  deleteColumns(numberColumns,delColumn);
+  delete [] delColumn;
+#else
+  defaultHandler_ = rhs->defaultHandler_;
+  if (defaultHandler_) 
+    handler_ = new CoinMessageHandler(*rhs->handler_);
+   else 
+    handler_ = rhs->handler_;
+  messages_ = rhs->messages_;
+  intParam_[ClpMaxNumIteration] = rhs->intParam_[ClpMaxNumIteration];
+  intParam_[ClpMaxNumIterationHotStart] = 
+    rhs->intParam_[ClpMaxNumIterationHotStart];
+
+  dblParam_[ClpDualObjectiveLimit] = rhs->dblParam_[ClpDualObjectiveLimit];
+  dblParam_[ClpPrimalObjectiveLimit] = rhs->dblParam_[ClpPrimalObjectiveLimit];
+  dblParam_[ClpDualTolerance] = rhs->dblParam_[ClpDualTolerance];
+  dblParam_[ClpPrimalTolerance] = rhs->dblParam_[ClpPrimalTolerance];
+  dblParam_[ClpObjOffset] = rhs->dblParam_[ClpObjOffset];
+  dblParam_[ClpMaxSeconds] = rhs->dblParam_[ClpMaxSeconds];
+
+  strParam_[ClpProbName] = rhs->strParam_[ClpProbName];
+
+  optimizationDirection_ = rhs->optimizationDirection_;
+  objectiveValue_=rhs->objectiveValue_;
+  smallElement_ = rhs->smallElement_;
+  numberIterations_ = rhs->numberIterations_;
+  solveType_ = rhs->solveType_;
+  problemStatus_ = rhs->problemStatus_;
+  // check valid lists
+  int numberBad=0;
+  int i;
+  for (i=0;i<numberRows;i++)
+    if (whichRow[i]<0||whichRow[i]>=rhs->numberRows_)
+      numberBad++;
+  if (numberBad)
+    throw CoinError("bad row list", "subproblem constructor", "ClpModel");
+  numberBad=0;
+  for (i=0;i<numberColumns;i++)
+    if (whichColumn[i]<0||whichColumn[i]>=rhs->numberColumns_)
+      numberBad++;
+  if (numberBad)
+    throw CoinError("bad column list", "subproblem constructor", "ClpModel");
+  numberRows_ = numberRows;
+  numberColumns_ = numberColumns;
+  if (!dropNames) {
+    unsigned int maxLength=0;
+    int iRow;
+    rowNames_.reserve(numberRows_);
+    for (iRow=0;iRow<numberRows_;iRow++) {
+      rowNames_[iRow] = rhs->rowNames_[whichRow[iRow]];
+      maxLength = max(maxLength,(unsigned int) strlen(rowNames_[iRow].c_str()));
+    }
+    int iColumn;
+    columnNames_.reserve(numberColumns_);
+    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
+      columnNames_[iColumn] = rhs->columnNames_[whichColumn[iColumn]];
+      maxLength = max(maxLength,(unsigned int) strlen(columnNames_[iColumn].c_str()));
+    }
+    lengthNames_=(int) maxLength;
+  } else {
+    lengthNames_ = 0;
+    rowNames_ = std::vector<std::string> ();
+    columnNames_ = std::vector<std::string> ();
+  }
+  if (rhs->integerType_&&!dropIntegers) {
+    integerType_ = whichChar(rhs->integerType_,numberColumns,whichColumn);
+  } else {
+    integerType_ = NULL;
+  }
+  if (rhs->rowActivity_) {
+    rowActivity_=whichDouble(rhs->rowActivity_,numberRows,whichRow);
+    dual_=whichDouble(rhs->dual_,numberRows,whichRow);
+    columnActivity_=whichDouble(rhs->columnActivity_,numberColumns,
+                                whichColumn);
+    reducedCost_=whichDouble(rhs->reducedCost_,numberColumns,
+                                whichColumn);
+  } else {
+    rowActivity_=NULL;
+    columnActivity_=NULL;
+    dual_=NULL;
+    reducedCost_=NULL;
+  }
+  rowLower_=whichDouble(rhs->rowLower_,numberRows,whichRow);
+  rowUpper_=whichDouble(rhs->rowUpper_,numberRows,whichRow);
+  columnLower_=whichDouble(rhs->columnLower_,numberColumns,whichColumn);
+  columnUpper_=whichDouble(rhs->columnUpper_,numberColumns,whichColumn);
+  if (rhs->objective_)
+    objective_  = rhs->objective_->subsetClone(numberColumns,whichColumn);
+  else
+    objective_ = NULL;
+  rowObjective_=whichDouble(rhs->rowObjective_,numberRows,whichRow);
+  // status has to be done in two stages
+  status_ = new unsigned char[numberColumns_+numberRows_];
+  unsigned char * rowStatus = whichUnsignedChar(rhs->status_+rhs->numberColumns_,
+                                                numberRows_,whichRow);
+  unsigned char * columnStatus = whichUnsignedChar(rhs->status_,
+                                                numberColumns_,whichColumn);
+  memcpy(status_+numberColumns_,rowStatus,numberRows_);
+  delete [] rowStatus;
+  memcpy(status_,columnStatus,numberColumns_);
+  delete [] columnStatus;
+  ray_ = NULL;
+  if (problemStatus_==1)
+    ray_ = whichDouble (rhs->ray_,numberRows,whichRow);
+  else if (problemStatus_==2)
+    ray_ = whichDouble (rhs->ray_,numberColumns,whichColumn);
+  if (rhs->rowCopy_) {
+    rowCopy_ = rhs->rowCopy_->subsetClone(numberRows,whichRow,
+                                          numberColumns,whichColumn);
+  } else {
+    rowCopy_=NULL;
+  }
+  if (rhs->quadraticObjective_) {
+    quadraticObjective_ = new CoinPackedMatrix(*rhs->quadraticObjective_,
+                                               numberColumns,whichColumn,
+                                               numberColumns,whichColumn);
+  } else {
+    quadraticObjective_=NULL;
+  }
+  matrix_=NULL;
+  if (rhs->matrix_) {
+    matrix_ = rhs->matrix_->subsetClone(numberRows,whichRow,
+                                        numberColumns,whichColumn);
+  }
+#endif
 }
 //#############################################################################

branches/pre/ClpNetworkBasis.cpp

@@ -561 +561 @@
 
 /* Updates one column (FTRAN) from region2 */
-int 
+double
 ClpNetworkBasis::updateColumn (  CoinIndexedVector * regionSparse,
-                                 CoinIndexedVector * regionSparse2)
+                                 CoinIndexedVector * regionSparse2,
+                                 int pivotRow)
 {
   regionSparse->clear (  );
@@ -573 +574 @@
   int i;
   bool doTwo = (numberNonZero==2);
-  int i0;
-  int i1;
+  int i0=-1;
+  int i1=-1;
   if (doTwo) {
     i0 = regionIndex2[0];
     i1 = regionIndex2[1];
-    doTwo =  (region2[i0]*region2[i1]<0.0);
-  }
-  if (doTwo) {
-    // If just +- 1 then could go backwards on depth until join
-    region[i0] = region2[i0];
-    region2[i0]=0.0;
-    region[i1] = region2[i1];
-    region2[i1]=0.0;
-    int iDepth0 = depth_[i0];
-    int iDepth1 = depth_[i1];
-    if (iDepth1>iDepth0) {
-      int temp = i0;
-      i0 = i1;
-      i1 = temp;
-      temp = iDepth0;
-      iDepth0 = iDepth1;
-      iDepth1 = temp;
-    }
-    numberNonZero=0;
-    while (iDepth0>iDepth1) {
-      double pivotValue = region[i0];
-      // put back now ?
-      int iBack = permuteBack_[i0];
-      regionIndex2[numberNonZero++] = iBack;
-      int otherRow = parent_[i0];
-      region2[iBack] = pivotValue*sign_[i0];
-      region[i0] =0.0;
-      region[otherRow] += pivotValue;
-      iDepth0--;
-      i0 = otherRow;
-    }
-    while (i0!=i1) {
-      double pivotValue = region[i0];
-      // put back now ?
-      int iBack = permuteBack_[i0];
-      regionIndex2[numberNonZero++] = iBack;
-      int otherRow = parent_[i0];
-      region2[iBack] = pivotValue*sign_[i0];
-      region[i0] =0.0;
-      region[otherRow] += pivotValue;
-      i0 = otherRow;
-      double pivotValue1 = region[i1];
-      // put back now ?
-      int iBack1 = permuteBack_[i1];
-      regionIndex2[numberNonZero++] = iBack1;
-      int otherRow1 = parent_[i1];
-      region2[iBack1] = pivotValue1*sign_[i1];
-      region[i1] =0.0;
-      region[otherRow1] += pivotValue1;
-      i1 = otherRow1;
-    }
-  } else {
-    // set up linked lists at each depth
-    // stack2 is start, stack is next
-    int greatestDepth=-1;
-    //mark_[numberRows_]=1;
-    for (i=0;i<numberNonZero;i++) {
-      int j = regionIndex2[i];
-      double value = region2[j];
-      region2[j]=0.0;
-      region[j]=value;
-      regionIndex[i]=j;
-      int iDepth = depth_[j];
-      if (iDepth>greatestDepth) 
-        greatestDepth = iDepth;
-      // and back until marked
-      while (!mark_[j]) {
-        int iNext = stack2_[iDepth];
-        stack2_[iDepth]=j;
-        stack_[j]=iNext;
-        mark_[j]=1;
-        iDepth--;
-        j=parent_[j];
-      }
-    }
-    numberNonZero=0;
-    for (;greatestDepth>=0; greatestDepth--) {
-      int iPivot = stack2_[greatestDepth];
-      stack2_[greatestDepth]=-1;
-      while (iPivot>=0) {
-        mark_[iPivot]=0;
-        double pivotValue = region[iPivot];
-        if (pivotValue) {
+  }
+  double returnValue=0.0;
+  bool packed = regionSparse2->packedMode();
+  if (packed) {
+    if (doTwo&&region2[0]*region2[1]<0.0) {
+      region[i0] = region2[0];
+      region2[0]=0.0;
+      region[i1] = region2[1];
+      region2[1]=0.0;
+      int iDepth0 = depth_[i0];
+      int iDepth1 = depth_[i1];
+      if (iDepth1>iDepth0) {
+        int temp = i0;
+        i0 = i1;
+        i1 = temp;
+        temp = iDepth0;
+        iDepth0 = iDepth1;
+        iDepth1 = temp;
+      }
+      numberNonZero=0;
+      if (pivotRow<0) {
+        while (iDepth0>iDepth1) {
+          double pivotValue = region[i0];
           // put back now ?
-          int iBack = permuteBack_[iPivot];
+          int iBack = permuteBack_[i0];
+          region2[numberNonZero] = pivotValue*sign_[i0];
           regionIndex2[numberNonZero++] = iBack;
-          int otherRow = parent_[iPivot];
-          region2[iBack] = pivotValue*sign_[iPivot];
-          region[iPivot] =0.0;
+          int otherRow = parent_[i0];
+          region[i0] =0.0;
           region[otherRow] += pivotValue;
-        }
-        iPivot = stack_[iPivot];
-      }
-    }
+          iDepth0--;
+          i0 = otherRow;
+        }
+        while (i0!=i1) {
+          double pivotValue = region[i0];
+          // put back now ?
+          int iBack = permuteBack_[i0];
+          region2[numberNonZero] = pivotValue*sign_[i0];
+          regionIndex2[numberNonZero++] = iBack;
+          int otherRow = parent_[i0];
+          region[i0] =0.0;
+          region[otherRow] += pivotValue;
+          i0 = otherRow;
+          double pivotValue1 = region[i1];
+          // put back now ?
+          int iBack1 = permuteBack_[i1];
+          region2[numberNonZero] = pivotValue1*sign_[i1];
+          regionIndex2[numberNonZero++] = iBack1;
+          int otherRow1 = parent_[i1];
+          region[i1] =0.0;
+          region[otherRow1] += pivotValue1;
+          i1 = otherRow1;
+        }
+      } else {
+        while (iDepth0>iDepth1) {
+          double pivotValue = region[i0];
+          // put back now ?
+          int iBack = permuteBack_[i0];
+          double value = pivotValue*sign_[i0];
+          region2[numberNonZero] = value;
+          regionIndex2[numberNonZero++] = iBack;
+          if (iBack==pivotRow)
+            returnValue = value;
+          int otherRow = parent_[i0];
+          region[i0] =0.0;
+          region[otherRow] += pivotValue;
+          iDepth0--;
+          i0 = otherRow;
+        }
+        while (i0!=i1) {
+          double pivotValue = region[i0];
+          // put back now ?
+          int iBack = permuteBack_[i0];
+          double value = pivotValue*sign_[i0];
+          region2[numberNonZero] = value;
+          regionIndex2[numberNonZero++] = iBack;
+          if (iBack==pivotRow)
+            returnValue = value;
+          int otherRow = parent_[i0];
+          region[i0] =0.0;
+          region[otherRow] += pivotValue;
+          i0 = otherRow;
+          double pivotValue1 = region[i1];
+          // put back now ?
+          int iBack1 = permuteBack_[i1];
+          value = pivotValue1*sign_[i1];
+          region2[numberNonZero] = value;
+          regionIndex2[numberNonZero++] = iBack1;
+          if (iBack1==pivotRow)
+            returnValue = value;
+          int otherRow1 = parent_[i1];
+          region[i1] =0.0;
+          region[otherRow1] += pivotValue1;
+          i1 = otherRow1;
+        }
+      }
+    } else {
+      // set up linked lists at each depth
+      // stack2 is start, stack is next
+      int greatestDepth=-1;
+      //mark_[numberRows_]=1;
+      for (i=0;i<numberNonZero;i++) {
+        int j = regionIndex2[i];
+        double value = region2[j];
+        region2[j]=0.0;
+        region[j]=value;
+        regionIndex[i]=j;
+        int iDepth = depth_[j];
+        if (iDepth>greatestDepth) 
+          greatestDepth = iDepth;
+        // and back until marked
+        while (!mark_[j]) {
+          int iNext = stack2_[iDepth];
+          stack2_[iDepth]=j;
+          stack_[j]=iNext;
+          mark_[j]=1;
+          iDepth--;
+          j=parent_[j];
+        }
+      }
+      numberNonZero=0;
+      if (pivotRow<0) {
+        for (;greatestDepth>=0; greatestDepth--) {
+          int iPivot = stack2_[greatestDepth];
+          stack2_[greatestDepth]=-1;
+          while (iPivot>=0) {
+            mark_[iPivot]=0;
+            double pivotValue = region[iPivot];
+            if (pivotValue) {
+              // put back now ?
+              int iBack = permuteBack_[iPivot];
+              region2[numberNonZero] = pivotValue*sign_[iPivot];
+              regionIndex2[numberNonZero++] = iBack;
+              int otherRow = parent_[iPivot];
+              region[iPivot] =0.0;
+            region[otherRow] += pivotValue;
+            }
+            iPivot = stack_[iPivot];
+          }
+        }
+      } else {
+        for (;greatestDepth>=0; greatestDepth--) {
+          int iPivot = stack2_[greatestDepth];
+          stack2_[greatestDepth]=-1;
+          while (iPivot>=0) {
+            mark_[iPivot]=0;
+            double pivotValue = region[iPivot];
+            if (pivotValue) {
+              // put back now ?
+              int iBack = permuteBack_[iPivot];
+              double value = pivotValue*sign_[iPivot];
+              region2[numberNonZero] = value;
+              regionIndex2[numberNonZero++] = iBack;
+              if (iBack==pivotRow)
+                returnValue = value;
+              int otherRow = parent_[iPivot];
+              region[iPivot] =0.0;
+              region[otherRow] += pivotValue;
+            }
+            iPivot = stack_[iPivot];
+          }
+        }
+      }
+    }
+  } else {
+    if (doTwo&&region2[i0]*region2[i1]<0.0) {
+      // If just +- 1 then could go backwards on depth until join
+      region[i0] = region2[i0];
+      region2[i0]=0.0;
+      region[i1] = region2[i1];
+      region2[i1]=0.0;
+      int iDepth0 = depth_[i0];
+      int iDepth1 = depth_[i1];
+      if (iDepth1>iDepth0) {
+        int temp = i0;
+        i0 = i1;
+        i1 = temp;
+        temp = iDepth0;
+        iDepth0 = iDepth1;
+        iDepth1 = temp;
+      }
+      numberNonZero=0;
+      while (iDepth0>iDepth1) {
+        double pivotValue = region[i0];
+        // put back now ?
+        int iBack = permuteBack_[i0];
+        regionIndex2[numberNonZero++] = iBack;
+        int otherRow = parent_[i0];
+        region2[iBack] = pivotValue*sign_[i0];
+        region[i0] =0.0;
+        region[otherRow] += pivotValue;
+        iDepth0--;
+        i0 = otherRow;
+      }
+      while (i0!=i1) {
+        double pivotValue = region[i0];
+        // put back now ?
+        int iBack = permuteBack_[i0];
+        regionIndex2[numberNonZero++] = iBack;
+        int otherRow = parent_[i0];
+        region2[iBack] = pivotValue*sign_[i0];
+        region[i0] =0.0;
+        region[otherRow] += pivotValue;
+        i0 = otherRow;
+        double pivotValue1 = region[i1];
+        // put back now ?
+        int iBack1 = permuteBack_[i1];
+        regionIndex2[numberNonZero++] = iBack1;
+        int otherRow1 = parent_[i1];
+        region2[iBack1] = pivotValue1*sign_[i1];
+        region[i1] =0.0;
+        region[otherRow1] += pivotValue1;
+        i1 = otherRow1;
+      }
+    } else {
+      // set up linked lists at each depth
+      // stack2 is start, stack is next
+      int greatestDepth=-1;
+      //mark_[numberRows_]=1;
+      for (i=0;i<numberNonZero;i++) {
+        int j = regionIndex2[i];
+        double value = region2[j];
+        region2[j]=0.0;
+        region[j]=value;
+        regionIndex[i]=j;
+        int iDepth = depth_[j];
+        if (iDepth>greatestDepth) 
+          greatestDepth = iDepth;
+        // and back until marked
+        while (!mark_[j]) {
+          int iNext = stack2_[iDepth];
+          stack2_[iDepth]=j;
+          stack_[j]=iNext;
+          mark_[j]=1;
+          iDepth--;
+          j=parent_[j];
+        }
+      }
+      numberNonZero=0;
+      for (;greatestDepth>=0; greatestDepth--) {
+        int iPivot = stack2_[greatestDepth];
+        stack2_[greatestDepth]=-1;
+        while (iPivot>=0) {
+          mark_[iPivot]=0;
+          double pivotValue = region[iPivot];
+          if (pivotValue) {
+            // put back now ?
+            int iBack = permuteBack_[iPivot];
+            regionIndex2[numberNonZero++] = iBack;
+            int otherRow = parent_[iPivot];
+            region2[iBack] = pivotValue*sign_[iPivot];
+            region[iPivot] =0.0;
+            region[otherRow] += pivotValue;
+          }
+          iPivot = stack_[iPivot];
+        }
+      }
+    }
+    if (pivotRow>=0)
+      returnValue = region2[pivotRow];
   }
   region[numberRows_]=0.0;
@@ -684 +851 @@
  }
 #endif
-  return numberNonZero;
+  return returnValue;
 }
 /* Updates one column (FTRAN) to/from array 
@@ -832 +999 @@
   int i;
   int numberNonZero=0;
-  for (i=0;i<numberNonZero2;i++) {
-    int k = regionIndex2[i];
-    int j = permute_[k];
-    double value = region2[k];
-    region2[k]=0.0;
-    region[j]=value;
-    mark_[j]=1;
-    regionIndex[numberNonZero++]=j;
-  }
-  // copy back
-  // set up linked lists at each depth
-  // stack2 is start, stack is next
-  int greatestDepth=-1;
-  int smallestDepth=numberRows_;
-  //mark_[numberRows_]=1;
-  for (i=0;i<numberNonZero2;i++) {
-    int j = regionIndex[i];
-    double value = region[j];
-    region[j]=0.0;
-    region2[j]=value;
-    regionIndex2[i]=j;
-    // add in
-    int iDepth = depth_[j];
-    smallestDepth = min(iDepth,smallestDepth) ;
-    greatestDepth = max(iDepth,greatestDepth) ;
-    int jNext = stack2_[iDepth];
-    stack2_[iDepth]=j;
-    stack_[j]=jNext;
-    // and put all descendants on list
-    int iChild = descendant_[j];
-    while (iChild>=0) {
-      if (!mark_[iChild]) {
-        regionIndex2[numberNonZero++] = iChild;
-        mark_[iChild]=1;
-      }
-      iChild = rightSibling_[iChild];
-    }
-  }
-  for (;i<numberNonZero;i++) {
-    int j = regionIndex2[i];
-    // add in
-    int iDepth = depth_[j];
-    smallestDepth = min(iDepth,smallestDepth) ;
-    greatestDepth = max(iDepth,greatestDepth) ;
-    int jNext = stack2_[iDepth];
-    stack2_[iDepth]=j;
-    stack_[j]=jNext;
-    // and put all descendants on list
-    int iChild = descendant_[j];
-    while (iChild>=0) {
-      if (!mark_[iChild]) {
-        regionIndex2[numberNonZero++] = iChild;
-        mark_[iChild]=1;
-      }
-      iChild = rightSibling_[iChild];
-    }
-  }
-  numberNonZero2=0;
-  region2[numberRows_]=0.0;
-  int iDepth;
-  for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
-    int iPivot = stack2_[iDepth];
-    stack2_[iDepth]=-1;
-    while (iPivot>=0) {
-      mark_[iPivot]=0;
-      double pivotValue = region2[iPivot];
-      int otherRow = parent_[iPivot];
-      double otherValue = region2[otherRow];
-      pivotValue = sign_[iPivot]*pivotValue+otherValue;
-      region2[iPivot]=pivotValue;
-      if (pivotValue) 
-        regionIndex2[numberNonZero2++]=iPivot;
-      iPivot = stack_[iPivot];
+  bool packed = regionSparse2->packedMode();
+  if (packed) {
+    for (i=0;i<numberNonZero2;i++) {
+      int k = regionIndex2[i];
+      int j = permute_[k];
+      double value = region2[i];
+      region2[i]=0.0;
+      region[j]=value;
+      mark_[j]=1;
+      regionIndex[numberNonZero++]=j;
+    }
+    // set up linked lists at each depth
+    // stack2 is start, stack is next
+    int greatestDepth=-1;
+    int smallestDepth=numberRows_;
+    //mark_[numberRows_]=1;
+    for (i=0;i<numberNonZero2;i++) {
+      int j = regionIndex[i];
+      regionIndex2[i]=j;
+      // add in
+      int iDepth = depth_[j];
+      smallestDepth = min(iDepth,smallestDepth) ;
+      greatestDepth = max(iDepth,greatestDepth) ;
+      int jNext = stack2_[iDepth];
+      stack2_[iDepth]=j;
+      stack_[j]=jNext;
+      // and put all descendants on list
+      int iChild = descendant_[j];
+      while (iChild>=0) {
+        if (!mark_[iChild]) {
+          regionIndex2[numberNonZero++] = iChild;
+          mark_[iChild]=1;
+        }
+        iChild = rightSibling_[iChild];
+      }
+    }
+    for (;i<numberNonZero;i++) {
+      int j = regionIndex2[i];
+      // add in
+      int iDepth = depth_[j];
+      smallestDepth = min(iDepth,smallestDepth) ;
+      greatestDepth = max(iDepth,greatestDepth) ;
+      int jNext = stack2_[iDepth];
+      stack2_[iDepth]=j;
+      stack_[j]=jNext;
+      // and put all descendants on list
+      int iChild = descendant_[j];
+      while (iChild>=0) {
+        if (!mark_[iChild]) {
+          regionIndex2[numberNonZero++] = iChild;
+          mark_[iChild]=1;
+        }
+        iChild = rightSibling_[iChild];
+      }
+    }
+    numberNonZero2=0;
+    region[numberRows_]=0.0;
+    int iDepth;
+    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
+      int iPivot = stack2_[iDepth];
+      stack2_[iDepth]=-1;
+      while (iPivot>=0) {
+        mark_[iPivot]=0;
+        double pivotValue = region[iPivot];
+        int otherRow = parent_[iPivot];
+        double otherValue = region[otherRow];
+        pivotValue = sign_[iPivot]*pivotValue+otherValue;
+        region[iPivot]=pivotValue;
+        if (pivotValue) {
+          region2[numberNonZero2]=pivotValue;
+          regionIndex2[numberNonZero2++]=iPivot;
+        }
+        iPivot = stack_[iPivot];
+      }
+    }
+    // zero out
+    for (i=0;i<numberNonZero2;i++) {
+      int k = regionIndex2[i];
+      region[k]=0.0;
+    }
+  } else {
+    for (i=0;i<numberNonZero2;i++) {
+      int k = regionIndex2[i];
+      int j = permute_[k];
+      double value = region2[k];
+      region2[k]=0.0;
+      region[j]=value;
+      mark_[j]=1;
+      regionIndex[numberNonZero++]=j;
+    }
+    // copy back
+    // set up linked lists at each depth
+    // stack2 is start, stack is next
+    int greatestDepth=-1;
+    int smallestDepth=numberRows_;
+    //mark_[numberRows_]=1;
+    for (i=0;i<numberNonZero2;i++) {
+      int j = regionIndex[i];
+      double value = region[j];
+      region[j]=0.0;
+      region2[j]=value;
+      regionIndex2[i]=j;
+      // add in
+      int iDepth = depth_[j];
+      smallestDepth = min(iDepth,smallestDepth) ;
+      greatestDepth = max(iDepth,greatestDepth) ;
+      int jNext = stack2_[iDepth];
+      stack2_[iDepth]=j;
+      stack_[j]=jNext;
+      // and put all descendants on list
+      int iChild = descendant_[j];
+      while (iChild>=0) {
+        if (!mark_[iChild]) {
+          regionIndex2[numberNonZero++] = iChild;
+          mark_[iChild]=1;
+        }
+        iChild = rightSibling_[iChild];
+      }
+    }
+    for (;i<numberNonZero;i++) {
+      int j = regionIndex2[i];
+      // add in
+      int iDepth = depth_[j];
+      smallestDepth = min(iDepth,smallestDepth) ;
+      greatestDepth = max(iDepth,greatestDepth) ;
+      int jNext = stack2_[iDepth];
+      stack2_[iDepth]=j;
+      stack_[j]=jNext;
+      // and put all descendants on list
+      int iChild = descendant_[j];
+      while (iChild>=0) {
+        if (!mark_[iChild]) {
+          regionIndex2[numberNonZero++] = iChild;
+          mark_[iChild]=1;
+        }
+        iChild = rightSibling_[iChild];
+      }
+    }
+    numberNonZero2=0;
+    region2[numberRows_]=0.0;
+    int iDepth;
+    for (iDepth=smallestDepth;iDepth<=greatestDepth; iDepth++) {
+      int iPivot = stack2_[iDepth];
+      stack2_[iDepth]=-1;
+      while (iPivot>=0) {
+        mark_[iPivot]=0;
+        double pivotValue = region2[iPivot];
+        int otherRow = parent_[iPivot];
+        double otherValue = region2[otherRow];
+        pivotValue = sign_[iPivot]*pivotValue+otherValue;
+        region2[iPivot]=pivotValue;
+        if (pivotValue) 
+          regionIndex2[numberNonZero2++]=iPivot;
+        iPivot = stack_[iPivot];
+      }
     }
   }

branches/pre/ClpNetworkMatrix.cpp

@@ -369 +369 @@
   ClpPlusMinusOneMatrix* rowCopy =
     dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
+  bool packed = rowArray->packedMode();
   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     // do by column
     int iColumn;
-    double * markVector = y->denseVector(); // probably empty
+    assert (!y->getNumElements());
     CoinBigIndex j=0;
-    if (trueNetwork_) {
-      for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
-        double value = markVector[iColumn];
-        markVector[iColumn]=0.0;
-        int iRowM = indices_[j];
-        int iRowP = indices_[j+1];
-        value -= scalar*pi[iRowM];
-        value += scalar*pi[iRowP];
-        if (fabs(value)>zeroTolerance) {
-          index[numberNonZero++]=iColumn;
-          array[iColumn]=value;
-        }
+    if (packed) {
+      // need to expand pi into y
+      assert(y->capacity()>=numberRows);
+      double * piOld = pi;
+      pi = y->denseVector();
+      const int * whichRow = rowArray->getIndices();
+      int i;
+      // modify pi so can collapse to one loop
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=scalar*piOld[i];
+      }
+      if (trueNetwork_) {
+        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
+          double value = 0.0;
+          int iRowM = indices_[j];
+          int iRowP = indices_[j+1];
+          value -= pi[iRowM];
+          value += pi[iRowP];
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      } else {
+        // skip negative rows
+        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
+          double value = 0.0;
+          int iRowM = indices_[j];
+          int iRowP = indices_[j+1];
+          if (iRowM>=0)
+            value -= pi[iRowM];
+          if (iRowP>=0)
+            value += pi[iRowP];
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=0.0;
       }
     } else {
-      // skip negative rows
-      for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
-        double value = markVector[iColumn];
-        markVector[iColumn]=0.0;
-        int iRowM = indices_[j];
-        int iRowP = indices_[j+1];
-        if (iRowM>=0)
+      if (trueNetwork_) {
+        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
+          double value = 0.0;
+          int iRowM = indices_[j];
+          int iRowP = indices_[j+1];
           value -= scalar*pi[iRowM];
-        if (iRowP>=0)
           value += scalar*pi[iRowP];
-        if (fabs(value)>zeroTolerance) {
-          index[numberNonZero++]=iColumn;
-          array[iColumn]=value;
+          if (fabs(value)>zeroTolerance) {
+            index[numberNonZero++]=iColumn;
+            array[iColumn]=value;
+          }
+        }
+      } else {
+        // skip negative rows
+        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
+          double value = 0.0;
+          int iRowM = indices_[j];
+          int iRowP = indices_[j+1];
+          if (iRowM>=0)
+            value -= scalar*pi[iRowM];
+          if (iRowP>=0)
+            value += scalar*pi[iRowP];
+          if (fabs(value)>zeroTolerance) {
+            index[numberNonZero++]=iColumn;
+            array[iColumn]=value;
+          }
         }
       }
     }
     columnArray->setNumElements(numberNonZero);
-    y->setNumElements(0);
   } else {
     // do by row
@@ -430 +474 @@
   int numberToDo = y->getNumElements();
   const int * which = y->getIndices();
-  if (trueNetwork_) {
-    for (jColumn=0;jColumn<numberToDo;jColumn++) {
-      int iColumn = which[jColumn];
-      double value = 0.0;
-      CoinBigIndex j=iColumn<<1;
-      int iRowM = indices_[j];
-      int iRowP = indices_[j+1];
-      value -= pi[iRowM];
-      value += pi[iRowP];
-      if (fabs(value)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn]=value;
-      }
+  bool packed = rowArray->packedMode();
+  if (packed) {
+    // Must line up with y
+    // need to expand pi into y
+    int numberInRowArray = rowArray->getNumElements();
+    assert(y->capacity()>=model->numberRows());
+    double * piOld = pi;
+    pi = y->denseVector();
+    const int * whichRow = rowArray->getIndices();
+    int i;
+    for (i=0;i<numberInRowArray;i++) {
+      int iRow = whichRow[i];
+      pi[iRow]=piOld[i];
+    }
+    if (trueNetwork_) {
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        value -= pi[iRowM];
+        value += pi[iRowP];
+        array[jColumn]=value;
+      }
+    } else {
+      // skip negative rows
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        if (iRowM>=0)
+          value -= pi[iRowM];
+        if (iRowP>=0)
+          value += pi[iRowP];
+        array[jColumn]=value;
+      }
+    }
+    for (i=0;i<numberInRowArray;i++) {
+      int iRow = whichRow[i];
+      pi[iRow]=0.0;
     }
   } else {
-    // skip negative rows
-    for (jColumn=0;jColumn<numberToDo;jColumn++) {
-      int iColumn = which[jColumn];
-      double value = 0.0;
-      CoinBigIndex j=iColumn<<1;
-      int iRowM = indices_[j];
-      int iRowP = indices_[j+1];
-      if (iRowM>=0)
+    if (trueNetwork_) {
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
         value -= pi[iRowM];
-      if (iRowP>=0)
         value += pi[iRowP];
-      if (fabs(value)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn]=value;
+        if (fabs(value)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn]=value;
+        }
+      }
+    } else {
+      // skip negative rows
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        if (iRowM>=0)
+          value -= pi[iRowM];
+        if (iRowP>=0)
+          value += pi[iRowP];
+        if (fabs(value)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn]=value;
+        }
       }
     }
@@ -542 +632 @@
   return numberElements;
 }
+/* If element NULL returns number of elements in column part of basis,
+   If not NULL fills in as well */
+CoinBigIndex 
+ClpNetworkMatrix::fillBasis(const ClpSimplex * model,
+                                 const int * whichColumn, 
+                                 int numberBasic,
+                                 int numberColumnBasic,
+                                 int * indexRowU, int * indexColumnU,
+                                 double * elementU) const 
+{
+  int i;
+  CoinBigIndex numberElements=0;
+  if (elementU!=NULL) {
+    if (trueNetwork_) {
+      for (i=0;i<numberColumnBasic;i++) {
+        int iColumn = whichColumn[i];
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        indexRowU[numberElements]=iRowM;
+        indexColumnU[numberElements]=numberBasic;
+        elementU[numberElements]=-1.0;
+        indexRowU[numberElements+1]=iRowP;
+        indexColumnU[numberElements+1]=numberBasic;
+        elementU[numberElements+1]=1.0;
+        numberElements+=2;
+        numberBasic++;
+      }
+    } else {
+      for (i=0;i<numberColumnBasic;i++) {
+        int iColumn = whichColumn[i];
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        if (iRowM>=0) {
+          indexRowU[numberElements]=iRowM;
+          indexColumnU[numberElements]=numberBasic;
+          elementU[numberElements++]=-1.0;
+        }
+        if (iRowP>=0) {
+          indexRowU[numberElements]=iRowP;
+          indexColumnU[numberElements]=numberBasic;
+          elementU[numberElements++]=1.0;
+        }
+        numberBasic++;
+      }
+    }
+  } else {
+    if (trueNetwork_) {
+      numberElements = 2*numberColumnBasic;
+    } else {
+      for (i=0;i<numberColumnBasic;i++) {
+        int iColumn = whichColumn[i];
+        CoinBigIndex j=iColumn<<1;
+        int iRowM = indices_[j];
+        int iRowP = indices_[j+1];
+        if (iRowM>=0)
+          numberElements ++;
+        if (iRowP>=0)
+          numberElements ++;
+      }
+    }
+  }
+  return numberElements;
+}
 /* Unpacks a column into an CoinIndexedvector
-      Note that model is NOT const.  Bounds and objective could
-      be modified if doing column generation */
+ */
 void 
 ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
@@ -557 +711 @@
     rowArray->add(iRowP,1.0);
 }
+/* Unpacks a column into an CoinIndexedvector
+** in packed foramt
+Note that model is NOT const.  Bounds and objective could
+be modified if doing column generation (just for this variable) */
+void 
+ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
+                            CoinIndexedVector * rowArray,
+                            int iColumn) const
+{
+  int * index = rowArray->getIndices();
+  double * array = rowArray->denseVector();
+  int number = 0;
+  CoinBigIndex j=iColumn<<1;
+  int iRowM = indices_[j];
+  int iRowP = indices_[j+1];
+  if (iRowM>=0) {
+    array[number]=-1.0;
+    index[number++]=iRowM;
+  }
+  if (iRowP>=0) {
+    array[number]=1.0;
+    index[number++]=iRowP;
+  }
+  rowArray->setNumElements(number);
+  rowArray->setPackedMode(true);
+}
 /* Adds multiple of a column into an CoinIndexedvector
       You can use quickAdd to add to vector */

branches/pre/ClpNonLinearCost.cpp

@@ -364 +364 @@
 // purpose which will come to me later
 void 
-ClpNonLinearCost::checkInfeasibilities(bool toNearest)
+ClpNonLinearCost::checkInfeasibilities(double oldTolerance)
 {
   numberInfeasibilities_=0;
@@ -372 +372 @@
   sumInfeasibilities_ = 0.0;
   double primalTolerance = model_->currentPrimalTolerance();
-  
   int iSequence;
   double * solution = model_->solutionRegion();
@@ -378 +377 @@
   double * lower = model_->lowerRegion();
   double * cost = model_->costRegion();
+  bool toNearest = oldTolerance<=0.0;
     
   // nonbasic should be at a valid bound
@@ -454 +454 @@
       if (!toNearest) {
         // With increasing tolerances - we may be at wrong place
-        if (fabs(value-upperValue)>primalTolerance*1.0001) {
-          assert(fabs(value-lowerValue)<=primalTolerance*1.0001); 
+        if (fabs(value-upperValue)>oldTolerance*1.0001) {
+          assert(fabs(value-lowerValue)<=oldTolerance*1.0001); 
+          if  (fabs(value-lowerValue)>primalTolerance)
+            solution[iSequence]=lowerValue;
           model_->setStatus(iSequence,ClpSimplex::atLowerBound);
+        } else if  (fabs(value-upperValue)>primalTolerance) {
+          solution[iSequence]=upperValue;
         }
       } else {
@@ -470 +474 @@
       if (!toNearest) {
         // With increasing tolerances - we may be at wrong place
-        if (fabs(value-lowerValue)>primalTolerance*1.0001) {
-          assert(fabs(value-upperValue)<=primalTolerance*1.0001); 
-          model_->setStatus(iSequence,ClpSimplex::atUpperBound);
+        if (fabs(value-lowerValue)>oldTolerance*1.0001) {
+          assert(fabs(value-upperValue)<=oldTolerance*1.0001); 
+          if  (fabs(value-upperValue)>primalTolerance)
+            solution[iSequence]=upperValue;
+          model_->setStatus(iSequence,ClpSimplex::atLowerBound);
+        } else if  (fabs(value-lowerValue)>primalTolerance) {
+          solution[iSequence]=lowerValue;
         }
       } else {

branches/pre/ClpObjective.cpp

@@ -47 +47 @@
   return *this;
 }
+/* Subset clone.  Duplicates are allowed
+   and order is as given.
+*/
+ClpObjective * 
+ClpObjective::subsetClone (int numberColumns, 
+                           const int * whichColumns) const
+{
+  std::cerr<<"subsetClone not supported - ClpObjective"<<std::endl;
+  abort();
+}
+

branches/pre/ClpPackedMatrix.cpp

@@ -90 +90 @@
 {
   return new ClpPackedMatrix(*this);
+}
+/* Subset clone (without gaps).  Duplicates are allowed
+   and order is as given */
+ClpMatrixBase * 
+ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
+                              int numberColumns, 
+                              const int * whichColumns) const 
+{
+  return new ClpPackedMatrix(*this, numberRows, whichRows,
+                                   numberColumns, whichColumns);
+}
+/* Subset constructor (without gaps).  Duplicates are allowed
+   and order is as given */
+ClpPackedMatrix::ClpPackedMatrix (
+                       const ClpPackedMatrix & rhs,
+                       int numberRows, const int * whichRows,
+                       int numberColumns, const int * whichColumns)
+: ClpMatrixBase(rhs)
+{
+  matrix_ = new CoinPackedMatrix(*(rhs.matrix_),numberRows,whichRows,
+                                 numberColumns,whichColumns);
+  zeroElements_ = rhs.zeroElements_;
+}
+ClpPackedMatrix::ClpPackedMatrix (
+                       const CoinPackedMatrix & rhs,
+                       int numberRows, const int * whichRows,
+                       int numberColumns, const int * whichColumns)
+: ClpMatrixBase()
+{
+  matrix_ = new CoinPackedMatrix(rhs,numberRows,whichRows,
+                                 numberColumns,whichColumns);
+  zeroElements_ = false;
+  setType(1);
 }
 
@@ -227 +260 @@
   ClpPackedMatrix* rowCopy =
     dynamic_cast< ClpPackedMatrix*>(model->rowCopy());
+  bool packed = rowArray->packedMode();
   if (numberInRowArray>0.333*numberRows||!rowCopy) {
     // do by column
@@ -238 +272 @@
     int numberColumns = model->numberColumns();
     if (!y->getNumElements()) {
-      if (!rowScale) {
-        if (scalar==-1.0) {
+      if (packed) {
+        // need to expand pi into y
+        assert(y->capacity()>=numberRows);
+        double * piOld = pi;
+        pi = y->denseVector();
+        const int * whichRow = rowArray->getIndices();
+        int i;
+        if (!rowScale) {
+          // modify pi so can collapse to one loop
+          for (i=0;i<numberInRowArray;i++) {
+            int iRow = whichRow[i];
+            pi[iRow]=scalar*piOld[i];
+          }
           for (iColumn=0;iColumn<numberColumns;iColumn++) {
             double value = 0.0;
@@ -249 +294 @@
             }
             if (fabs(value)>zeroTolerance) {
+              array[numberNonZero]=value;
               index[numberNonZero++]=iColumn;
-              array[iColumn]=-value;
-            }
-          }
-        } else if (scalar==1.0) {
-          for (iColumn=0;iColumn<numberColumns;iColumn++) {
-            double value = 0.0;
-            CoinBigIndex j;
-            for (j=columnStart[iColumn];
-                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
-              int iRow = row[j];
-              value += pi[iRow]*elementByColumn[j];
-            }
-            if (fabs(value)>zeroTolerance) {
-              index[numberNonZero++]=iColumn;
-              array[iColumn]=value;
             }
           }
         } else {
-          for (iColumn=0;iColumn<numberColumns;iColumn++) {
-            double value = 0.0;
-            CoinBigIndex j;
-            for (j=columnStart[iColumn];
-                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
-              int iRow = row[j];
-              value += pi[iRow]*elementByColumn[j];
-            }
-            value *= scalar;
-            if (fabs(value)>zeroTolerance) {
-              index[numberNonZero++]=iColumn;
-              array[iColumn]=value;
-            }
-          }
-        }
-      } else {
-        // scaled
-        if (scalar==-1.0) {
+          // scaled
+          // modify pi so can collapse to one loop
+          for (i=0;i<numberInRowArray;i++) {
+            int iRow = whichRow[i];
+            pi[iRow]=scalar*piOld[i]*rowScale[iRow];
+          }
           for (iColumn=0;iColumn<numberColumns;iColumn++) {
             double value = 0.0;
@@ -293 +312 @@
                  j<columnStart[iColumn]+columnLength[iColumn];j++) {
               int iRow = row[j];
-              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+              value += pi[iRow]*elementByColumn[j];
             }
             value *= columnScale[iColumn];
             if (fabs(value)>zeroTolerance) {
+              array[numberNonZero]=value;
               index[numberNonZero++]=iColumn;
-              array[iColumn]=-value;
-            }
-          }
-        } else if (scalar==1.0) {
-          for (iColumn=0;iColumn<numberColumns;iColumn++) {
-            double value = 0.0;
-            CoinBigIndex j;
-            const double * columnScale = model->columnScale();
-            for (j=columnStart[iColumn];
-                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
-              int iRow = row[j];
-              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
-            }
-            value *= columnScale[iColumn];
-            if (fabs(value)>zeroTolerance) {
-              index[numberNonZero++]=iColumn;
-              array[iColumn]=value;
+            }
+          }
+        }
+        // zero out
+        for (i=0;i<numberInRowArray;i++) {
+          int iRow = whichRow[i];
+          pi[iRow]=0.0;
+        }
+      } else {
+        if (!rowScale) {
+          if (scalar==-1.0) {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j];
+              }
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=-value;
+              }
+            }
+          } else if (scalar==1.0) {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j];
+              }
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=value;
+              }
+            }
+          } else {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j];
+              }
+              value *= scalar;
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=value;
+              }
             }
           }
         } else {
-          for (iColumn=0;iColumn<numberColumns;iColumn++) {
-            double value = 0.0;
-            CoinBigIndex j;
-            const double * columnScale = model->columnScale();
-            for (j=columnStart[iColumn];
-                 j<columnStart[iColumn]+columnLength[iColumn];j++) {
-              int iRow = row[j];
-              value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
-            }
-            value *= scalar*columnScale[iColumn];
-            if (fabs(value)>zeroTolerance) {
-              index[numberNonZero++]=iColumn;
-              array[iColumn]=value;
+          // scaled
+          if (scalar==-1.0) {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              const double * columnScale = model->columnScale();
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+              }
+              value *= columnScale[iColumn];
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=-value;
+              }
+            }
+          } else if (scalar==1.0) {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              const double * columnScale = model->columnScale();
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+              }
+              value *= columnScale[iColumn];
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=value;
+              }
+            }
+          } else {
+            for (iColumn=0;iColumn<numberColumns;iColumn++) {
+              double value = 0.0;
+              CoinBigIndex j;
+              const double * columnScale = model->columnScale();
+              for (j=columnStart[iColumn];
+                   j<columnStart[iColumn]+columnLength[iColumn];j++) {
+                int iRow = row[j];
+                value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+              }
+              value *= scalar*columnScale[iColumn];
+              if (fabs(value)>zeroTolerance) {
+                index[numberNonZero++]=iColumn;
+                array[iColumn]=value;
+              }
             }
           }
@@ -336 +426 @@
       }
     } else {
+      assert(!packed);
       double * markVector = y->denseVector(); // not empty
       if (!rowScale) {
@@ -381 +472 @@
     rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
   }
+  if (packed)
+    columnArray->setPackedMode(true);
+  if (0) {
+    columnArray->checkClean();
+    int numberNonZero=columnArray->getNumElements();;
+    int * index = columnArray->getIndices();
+    double * array = columnArray->denseVector();
+    int i;
+    for (i=0;i<numberNonZero;i++) {
+      int j=index[i];
+      double value;
+      if (packed)
+        value=array[i];
+      else
+        value=array[j];
+      printf("Ti %d %d %g\n",i,j,value);
+    }
+  }
 }
 /* Return <code>x * A + y</code> in <code>z</code>. 
@@ -402 +511 @@
   const double * element = getElements();
   const int * whichRow = rowArray->getIndices();
+  bool packed = rowArray->packedMode();
   if (numberInRowArray>2||y->getNumElements()) {
     // do by rows
     // ** Row copy is already scaled
     int iRow;
-    double * markVector = y->denseVector(); // probably empty .. but
     int * mark = y->getIndices();
     int numberOriginal=y->getNumElements();
     int i;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = mark[i];
-      index[i]=iColumn;
-      array[iColumn]=markVector[iColumn];
-      markVector[iColumn]=0.0;
-    }
-    numberNonZero=numberOriginal;
-    // and set up mark as char array
-    char * marked = (char *) markVector;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      marked[iColumn]=0;
-    }
-
-    for (i=0;i<numberInRowArray;i++) {
-      iRow = whichRow[i]; 
-      double value = pi[iRow]*scalar;
-      CoinBigIndex j;
+    if (packed) {
+      assert(!numberOriginal);
+      numberNonZero=0;
+      // and set up mark as char array
+      char * marked = (char *) (index+columnArray->capacity());
+      double * array2 = y->denseVector();
+#ifdef CLP_DEBUG
+      int numberColumns = model->numberColumns();
+      for (i=0;i<numberColumns;i++) {
+        assert(!marked[i]);
+        assert(!array2[i]);
+      }
+#endif      
+      for (i=0;i<numberInRowArray;i++) {
+        iRow = whichRow[i]; 
+        double value = pi[i]*scalar;
+        CoinBigIndex j;
+        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array2[iColumn] += value*element[j];
+        }
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        if (marked[iColumn]) {
+          double value = array2[iColumn];
+          array2[iColumn]=0.0;
+          marked[iColumn]=0;
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+    } else {
+      double * markVector = y->denseVector(); // probably empty .. but
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = mark[i];
+        index[i]=iColumn;
+        array[iColumn]=markVector[iColumn];
+        markVector[iColumn]=0.0;
+      }
+      numberNonZero=numberOriginal;
+      // and set up mark as char array
+      char * marked = (char *) markVector;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        marked[iColumn]=0;
+      }
+      
+      for (i=0;i<numberInRowArray;i++) {
+        iRow = whichRow[i]; 
+        double value = pi[iRow]*scalar;
+        CoinBigIndex j;
+        for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array[iColumn] += value*element[j];
+        }
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        marked[iColumn]=0;
+        if (fabs(array[iColumn])>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+        } else {
+          array[iColumn]=0.0;
+        }
+      }
+    }
+  } else if (numberInRowArray==2) {
+    // do by rows when two rows
+    int numberOriginal;
+    int i;
+    CoinBigIndex j;
+    numberNonZero=0;
+
+    double value;
+    if (packed) {
+      int iRow0 = whichRow[0]; 
+      int iRow1 = whichRow[1]; 
+      double pi0 = pi[0];
+      double pi1 = pi[1];
+      if (rowStart[iRow0+1]-rowStart[iRow0]>
+          rowStart[iRow1+1]-rowStart[iRow1]) {
+        // do one with fewer first
+        iRow0=iRow1;
+        iRow1=whichRow[0];
+        pi0=pi1;
+        pi1=pi[0];
+      }
+      // and set up mark as char array
+      char * marked = (char *) (index+columnArray->capacity());
+      int * lookup = y->getIndices();
+      value = pi0*scalar;
+      for (j=rowStart[iRow0];j<rowStart[iRow0+1];j++) {
+        int iColumn = column[j];
+        double value2 = value*element[j];
+        array[numberNonZero] = value2;
+        marked[iColumn]=1;
+        lookup[iColumn]=numberNonZero;
+        index[numberNonZero++]=iColumn;
+      }
+      numberOriginal = numberNonZero;
+      value = pi1*scalar;
+      for (j=rowStart[iRow1];j<rowStart[iRow1+1];j++) {
+        int iColumn = column[j];
+        double value2 = value*element[j];
+        // I am assuming no zeros in matrix
+        if (marked[iColumn]) {
+          int iLookup = lookup[iColumn];
+          array[iLookup] += value2;
+        } else {
+          if (fabs(value2)>zeroTolerance) {
+            array[numberNonZero] = value2;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+      // get rid of tiny values and zero out marked
+      int nDelete=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        marked[iColumn]=0;
+        if (fabs(array[i])<=zeroTolerance) 
+          nDelete++;
+      }
+      if (nDelete) {
+        numberOriginal=numberNonZero;
+        numberNonZero=0;
+        for (i=0;i<numberOriginal;i++) {
+          int iColumn = index[i];
+          double value = array[i];
+          array[i]=0.0;
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+    } else {
+      int iRow = whichRow[0]; 
+      value = pi[iRow]*scalar;
       for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
         int iColumn = column[j];
-        if (!marked[iColumn]) {
-          marked[iColumn]=1;
+        double value2 = value*element[j];
+        index[numberNonZero++]=iColumn;
+        array[iColumn] = value2;
+      }
+      iRow = whichRow[1]; 
+      value = pi[iRow]*scalar;
+      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
+        int iColumn = column[j];
+        double value2 = value*element[j];
+        // I am assuming no zeros in matrix
+        if (array[iColumn])
+          value2 += array[iColumn];
+        else
           index[numberNonZero++]=iColumn;
-        }
-        array[iColumn] += value*element[j];
-      }
-    }
-    // get rid of tiny values and zero out marked
-    numberOriginal=numberNonZero;
-    numberNonZero=0;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      marked[iColumn]=0;
-      if (fabs(array[iColumn])>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-      } else {
-        array[iColumn]=0.0;
-      }
-    }
-  } else if (numberInRowArray==2) {
-    // do by rows when two rows
-    int iRow;
-    int numberOriginal;
-    int i;
-    numberNonZero=0;
-
-    double value;
-    iRow = whichRow[0]; 
-    value = pi[iRow]*scalar;
-    CoinBigIndex j;
-    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
-      int iColumn = column[j];
-      double value2 = value*element[j];
-      index[numberNonZero++]=iColumn;
-      array[iColumn] = value2;
-    }
-    iRow = whichRow[1]; 
-    value = pi[iRow]*scalar;
-    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
-      int iColumn = column[j];
-      double value2 = value*element[j];
-      // I am assuming no zeros in matrix
-      if (array[iColumn])
-        value2 += array[iColumn];
-      else
-        index[numberNonZero++]=iColumn;
-      array[iColumn] = value2;
-    }
-    // get rid of tiny values and zero out marked
-    numberOriginal=numberNonZero;
-    numberNonZero=0;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      if (fabs(array[iColumn])>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-      } else {
-        array[iColumn]=0.0;
+        array[iColumn] = value2;
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        if (fabs(array[iColumn])>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+        } else {
+          array[iColumn]=0.0;
+        }
       }
     }
@@ -493 +709 @@
     int iRow=rowArray->getIndices()[0];
     numberNonZero=0;
-    double value = pi[iRow]*scalar;
     CoinBigIndex j;
-    for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
-      int iColumn = column[j];
-      double value2 = value*element[j];
-      if (fabs(value2)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn] = value2;
+    if (packed) {
+      double value = pi[0]*scalar;
+      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
+        int iColumn = column[j];
+        double value2 = value*element[j];
+        if (fabs(value2)>zeroTolerance) {
+          array[numberNonZero] = value2;
+          index[numberNonZero++]=iColumn;
+        }
+      }
+    } else {
+      double value = pi[iRow]*scalar;
+      for (j=rowStart[iRow];j<rowStart[iRow+1];j++) {
+        int iColumn = column[j];
+        double value2 = value*element[j];
+        if (fabs(value2)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn] = value2;
+        }
       }
     }
@@ -532 +760 @@
   int numberToDo = y->getNumElements();
   const int * which = y->getIndices();
-  if (!rowScale) {
-    for (jColumn=0;jColumn<numberToDo;jColumn++) {
-      int iColumn = which[jColumn];
-      double value = 0.0;
-      CoinBigIndex j;
-      for (j=columnStart[iColumn];
-           j<columnStart[iColumn]+columnLength[iColumn];j++) {
-        int iRow = row[j];
-        value += pi[iRow]*elementByColumn[j];
-      }
-      if (fabs(value)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn]=value;
-      }
+  bool packed = rowArray->packedMode();
+  if (packed) {
+    // need to expand pi into y
+    int numberInRowArray = rowArray->getNumElements();
+    assert(y->capacity()>=model->numberRows());
+    double * piOld = pi;
+    pi = y->denseVector();
+    const int * whichRow = rowArray->getIndices();
+    int i;
+    // Do NOT squash small elements - must line up with  y
+    if (!rowScale) {
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=piOld[i];
+      }
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j;
+        for (j=columnStart[iColumn];
+             j<columnStart[iColumn]+columnLength[iColumn];j++) {
+          int iRow = row[j];
+          value += pi[iRow]*elementByColumn[j];
+        }
+        array[jColumn]=value;
+      }
+    } else {
+      // scaled
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=rowScale[iRow]*piOld[i];
+      }
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j;
+        const double * columnScale = model->columnScale();
+        for (j=columnStart[iColumn];
+             j<columnStart[iColumn]+columnLength[iColumn];j++) {
+          int iRow = row[j];
+          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+        }
+        value *= columnScale[iColumn];
+        array[jColumn]=value;
+      }
+    }
+    for (i=0;i<numberInRowArray;i++) {
+      int iRow = whichRow[i];
+      pi[iRow]=0.0;
     }
   } else {
-    // scaled
-    for (jColumn=0;jColumn<numberToDo;jColumn++) {
-      int iColumn = which[jColumn];
-      double value = 0.0;
-      CoinBigIndex j;
-      const double * columnScale = model->columnScale();
-      for (j=columnStart[iColumn];
-           j<columnStart[iColumn]+columnLength[iColumn];j++) {
-        int iRow = row[j];
-        value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
-      }
-      value *= columnScale[iColumn];
-      if (fabs(value)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn]=value;
+    if (!rowScale) {
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j;
+        for (j=columnStart[iColumn];
+             j<columnStart[iColumn]+columnLength[iColumn];j++) {
+          int iRow = row[j];
+          value += pi[iRow]*elementByColumn[j];
+        }
+        if (fabs(value)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn]=value;
+        }
+      }
+    } else {
+      // scaled
+      for (jColumn=0;jColumn<numberToDo;jColumn++) {
+        int iColumn = which[jColumn];
+        double value = 0.0;
+        CoinBigIndex j;
+        const double * columnScale = model->columnScale();
+        for (j=columnStart[iColumn];
+             j<columnStart[iColumn]+columnLength[iColumn];j++) {
+          int iRow = row[j];
+          value += pi[iRow]*elementByColumn[j]*rowScale[iRow];
+        }
+        value *= columnScale[iColumn];
+        if (fabs(value)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn]=value;
+        }
       }
     }
@@ -681 +961 @@
         }
       }
+    }
+  }
+  return numberElements;
+}
+/* If element NULL returns number of elements in column part of basis,
+   If not NULL fills in as well */
+CoinBigIndex 
+ClpPackedMatrix::fillBasis(const ClpSimplex * model,
+                           const int * whichColumn, 
+                           int numberBasic,
+                           int numberColumnBasic,
+                           int * indexRowU, int * indexColumnU,
+                           double * elementU) const 
+{
+  const int * columnLength = matrix_->getVectorLengths(); 
+  int i;
+  CoinBigIndex numberElements=0;
+  if (elementU!=NULL) {
+    // fill
+    const CoinBigIndex * columnStart = matrix_->getVectorStarts();
+    const double * rowScale = model->rowScale();
+    const int * row = matrix_->getIndices();
+    const double * elementByColumn = matrix_->getElements();
+    if (!zeroElements_) {
+      if (!rowScale) {
+        // no scaling
+        for (i=0;i<numberColumnBasic;i++) {
+          int iColumn = whichColumn[i];
+          CoinBigIndex j;
+          for (j=columnStart[iColumn];
+               j<columnStart[iColumn]+columnLength[iColumn];j++) {
+            indexRowU[numberElements]=row[j];
+            indexColumnU[numberElements]=numberBasic;
+            elementU[numberElements++]=elementByColumn[j];
+          }
+          numberBasic++;
+        }
+      } else {
+        // scaling
+        const double * columnScale = model->columnScale();
+        for (i=0;i<numberColumnBasic;i++) {
+          int iColumn = whichColumn[i];
+          CoinBigIndex j;
+          double scale = columnScale[iColumn];
+          for (j=columnStart[iColumn];
+               j<columnStart[iColumn]+columnLength[iColumn];j++) {
+            int iRow = row[j];
+            indexRowU[numberElements]=iRow;
+            indexColumnU[numberElements]=numberBasic;
+            elementU[numberElements++]=
+              elementByColumn[j]*scale*rowScale[iRow];
+          }
+          numberBasic++;
+        }
+      }
+    } else {
+      // there are zero elements so need to look more closely
+      if (!rowScale) {
+        // no scaling
+        for (i=0;i<numberColumnBasic;i++) {
+          int iColumn = whichColumn[i];
+          CoinBigIndex j;
+          for (j=columnStart[iColumn];
+               j<columnStart[iColumn]+columnLength[iColumn];j++) {
+            double value = elementByColumn[j];
+            if (value) {
+              indexRowU[numberElements]=row[j];
+              indexColumnU[numberElements]=numberBasic;
+              elementU[numberElements++]=value;
+            }
+          }
+          numberBasic++;
+        }
+      } else {
+        // scaling
+        const double * columnScale = model->columnScale();
+        for (i=0;i<numberColumnBasic;i++) {
+          int iColumn = whichColumn[i];
+          CoinBigIndex j;
+          double scale = columnScale[iColumn];
+          for (j=columnStart[iColumn];
+               j<columnStart[iColumn]+columnLength[i];j++) {
+            double value = elementByColumn[j];
+            if (value) {
+              int iRow = row[j];
+              indexRowU[numberElements]=iRow;
+              indexColumnU[numberElements]=numberBasic;
+              elementU[numberElements++]=value*scale*rowScale[iRow];
+            }
+          }
+          numberBasic++;
+        }
+      }
+    }
+  } else {
+    // just count - can be over so ignore zero problem
+    for (i=0;i<numberColumnBasic;i++) {
+      int iColumn = whichColumn[i];
+      numberElements += columnLength[iColumn];
     }
   }
@@ -893 +1272 @@
 }
 /* Unpacks a column into an CoinIndexedvector
-      Note that model is NOT const.  Bounds and objective could
-      be modified if doing column generation */
+ */
 void 
 ClpPackedMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
@@ -918 +1296 @@
       rowArray->add(iRow,elementByColumn[i]*scale*rowScale[iRow]);
     }
+  }
+}
+/* Unpacks a column into an CoinIndexedvector
+** in packed foramt
+Note that model is NOT const.  Bounds and objective could
+be modified if doing column generation (just for this variable) */
+void 
+ClpPackedMatrix::unpackPacked(ClpSimplex * model,
+                            CoinIndexedVector * rowArray,
+                            int iColumn) const
+{
+  const double * rowScale = model->rowScale();
+  const int * row = matrix_->getIndices();
+  const CoinBigIndex * columnStart = matrix_->getVectorStarts();
+  const int * columnLength = matrix_->getVectorLengths(); 
+  const double * elementByColumn = matrix_->getElements();
+  CoinBigIndex i;
+  if (!rowScale) {
+    int j=columnStart[iColumn];
+    rowArray->createPacked(columnLength[iColumn],
+                           row+j,elementByColumn+j);
+  } else {
+    // apply scaling
+    double scale = model->columnScale()[iColumn];
+    int * index = rowArray->getIndices();
+    double * array = rowArray->denseVector();
+    int number = 0;
+    for (i=columnStart[iColumn];
+         i<columnStart[iColumn]+columnLength[iColumn];i++) {
+      int iRow = row[i];
+      array[number]=elementByColumn[i]*scale*rowScale[iRow];
+      index[number++]=iRow;
+    }
+    rowArray->setNumElements(number);
+    rowArray->setPackedMode(true);
   }
 }

branches/pre/ClpPlusMinusOneMatrix.cpp

@@ -51 +51 @@
   columnOrdered_=rhs.columnOrdered_;
   if (numberColumns_) {
-    indices_ = new int [ 2*numberColumns_];
-    memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
+    int numberElements = rhs.startPositive_[numberColumns_];
+    indices_ = new int [ numberElements];
+    memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
     startPositive_ = new int [ numberColumns_+1];
     memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
@@ -165 +166 @@
     columnOrdered_=rhs.columnOrdered_;
     if (numberColumns_) {
-      indices_ = new int [ 2*numberColumns_];
-      memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
+      int numberElements = rhs.startPositive_[numberColumns_];
+      indices_ = new int [ numberElements];
+      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
       startPositive_ = new int [ numberColumns_+1];
       memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
@@ -182 +184 @@
   return new ClpPlusMinusOneMatrix(*this);
 }
+/* Subset clone (without gaps).  Duplicates are allowed
+   and order is as given */
+ClpMatrixBase * 
+ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
+                              int numberColumns, 
+                              const int * whichColumns) const 
+{
+  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
+                                   numberColumns, whichColumns);
+}
+/* Subset constructor (without gaps).  Duplicates are allowed
+   and order is as given */
+ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
+                       const ClpPlusMinusOneMatrix & rhs,
+                       int numberRows, const int * whichRow,
+                       int numberColumns, const int * whichColumn)
+: ClpMatrixBase(rhs)
+{  
+  elements_ = NULL;
+  startPositive_ = NULL;
+  startNegative_ = NULL;
+  lengths_=NULL;
+  indices_=NULL;
+  numberRows_=0;
+  numberColumns_=0;
+  columnOrdered_=rhs.columnOrdered_;
+  if (numberRows<=0||numberColumns<=0) {
+    startPositive_ = new int[1];
+    startPositive_[0] = 0;
+  } else {
+    numberColumns_ = numberColumns;
+    numberRows_ = numberRows;
+    const int * index1 = rhs.indices_;
+    int * startPositive1 = rhs.startPositive_;
+
+    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
+    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
+    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
+    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
+    // Throw exception if rhs empty
+    if (numberMajor1 <= 0 || numberMinor1 <= 0)
+      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
+    // Array to say if an old row is in new copy
+    int * newRow = new int [numberMinor1];
+    int iRow;
+    for (iRow=0;iRow<numberMinor1;iRow++) 
+      newRow[iRow] = -1;
+    // and array for duplicating rows
+    int * duplicateRow = new int [numberMinor];
+    int numberBad=0;
+    for (iRow=0;iRow<numberMinor;iRow++) {
+      duplicateRow[iRow] = -1;
+      int kRow = whichRow[iRow];
+      if (kRow>=0  && kRow < numberMinor1) {
+        if (newRow[kRow]<0) {
+          // first time
+          newRow[kRow]=iRow;
+        } else {
+          // duplicate
+          int lastRow = newRow[kRow];
+          newRow[kRow]=iRow;
+          duplicateRow[iRow] = lastRow;
+        }
+      } else {
+        // bad row
+        numberBad++;
+      }
+    }
+
+    if (numberBad)
+      throw CoinError("bad minor entries", 
+                      "subset constructor", "ClpPlusMinusOneMatrix");
+    // now get size and check columns
+    int size = 0;
+    int iColumn;
+    numberBad=0;
+    for (iColumn=0;iColumn<numberMajor;iColumn++) {
+      int kColumn = whichColumn[iColumn];
+      if (kColumn>=0  && kColumn <numberMajor1) {
+        int i;
+        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
+          int kRow = index1[i];
+          kRow = newRow[kRow];
+          while (kRow>=0) {
+            size++;
+            kRow = duplicateRow[kRow];
+          }
+        }
+      } else {
+        // bad column
+        numberBad++;
+      }
+    }
+    if (numberBad)
+      throw CoinError("bad major entries", 
+                      "subset constructor", "ClpPlusMinusOneMatrix");
+    // now create arrays
+    startPositive_ = new int [numberMajor+1];
+    startNegative_ = new int [numberMajor];
+    indices_ = new int[size];
+    // and fill them
+    size = 0;
+    startPositive_[0]=0;
+    int * startNegative1 = rhs.startNegative_;
+    for (iColumn=0;iColumn<numberMajor;iColumn++) {
+      int kColumn = whichColumn[iColumn];
+      int i;
+      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
+        int kRow = index1[i];
+        kRow = newRow[kRow];
+        while (kRow>=0) {
+          indices_[size++] = kRow;
+          kRow = duplicateRow[kRow];
+        }
+      }
+      startNegative_[iColumn] = size;
+      for (;i<startPositive1[kColumn+1];i++) {
+        int kRow = index1[i];
+        kRow = newRow[kRow];
+        while (kRow>=0) {
+          indices_[size++] = kRow;
+          kRow = duplicateRow[kRow];
+        }
+      }
+      startPositive_[iColumn+1] = size;
+    }
+    delete [] newRow;
+    delete [] duplicateRow;
+  }
+}
+
 
 /* Returns a new matrix in reverse order without gaps */
@@ -323 +456 @@
   double zeroTolerance = model->factorization()->zeroTolerance();
   int numberRows = model->numberRows();
+  bool packed = rowArray->packedMode();
   ClpPlusMinusOneMatrix* rowCopy =
     dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
   if (numberInRowArray>0.3*numberRows||!rowCopy) {
+    assert (!y->getNumElements());
     // do by column
+    // Need to expand if packed mode
     int iColumn;
-    double * markVector = y->denseVector(); // probably empty
     CoinBigIndex j=0;
     assert (columnOrdered_);
-    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
-      double value2 = 0.0;
-      for (;j<startNegative_[iColumn];j++) {
-        int iRow = indices_[j];
-        value2 += pi[iRow];
-      }
-      for (;j<startPositive_[iColumn+1];j++) {
-        int iRow = indices_[j];
-        value2 -= pi[iRow];
-      }
-      double value = markVector[iColumn];
-      markVector[iColumn]=0.0;
-      value += scalar*value2;
-      if (fabs(value)>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-        array[iColumn]=value;
+    if (packed) {
+      // need to expand pi into y
+      assert(y->capacity()>=numberRows);
+      double * piOld = pi;
+      pi = y->denseVector();
+      const int * whichRow = rowArray->getIndices();
+      int i;
+      // modify pi so can collapse to one loop
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=scalar*piOld[i];
+      }
+      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
+        double value = 0.0;
+        for (;j<startNegative_[iColumn];j++) {
+          int iRow = indices_[j];
+          value += pi[iRow];
+        }
+        for (;j<startPositive_[iColumn+1];j++) {
+          int iRow = indices_[j];
+          value -= pi[iRow];
+        }
+        if (fabs(value)>zeroTolerance) {
+          array[numberNonZero]=value;
+          index[numberNonZero++]=iColumn;
+        }
+      }
+      for (i=0;i<numberInRowArray;i++) {
+        int iRow = whichRow[i];
+        pi[iRow]=0.0;
+      }
+    } else {
+      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
+        double value = 0.0;
+        for (;j<startNegative_[iColumn];j++) {
+          int iRow = indices_[j];
+          value += pi[iRow];
+        }
+        for (;j<startPositive_[iColumn+1];j++) {
+          int iRow = indices_[j];
+          value -= pi[iRow];
+        }
+        value += scalar;
+        if (fabs(value)>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+          array[iColumn]=value;
+        }
       }
     }
     columnArray->setNumElements(numberNonZero);
-    y->setNumElements(0);
   } else {
     // do by row
@@ -376 +541 @@
   const CoinBigIndex * startNegative = startNegative_;
   const int * whichRow = rowArray->getIndices();
+  bool packed = rowArray->packedMode();
   if (numberInRowArray>2||y->getNumElements()) {
     // do by rows
@@ -383 +549 @@
     int numberOriginal=y->getNumElements();
     int i;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = mark[i];
-      index[i]=iColumn;
-      array[iColumn]=markVector[iColumn];
-      markVector[iColumn]=0.0;
-    }
-    numberNonZero=numberOriginal;
-    // and set up mark as char array
-    char * marked = (char *) markVector;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      marked[iColumn]=0;
-    }
-    for (i=0;i<numberInRowArray;i++) {
-      iRow = whichRow[i]; 
-      double value = pi[iRow]*scalar;
-      CoinBigIndex j;
-      for (j=startPositive[iRow];j<startNegative[iRow];j++) {
-        int iColumn = column[j];
-        if (!marked[iColumn]) {
-          marked[iColumn]=1;
+    if (packed) {
+      assert(!numberOriginal);
+      numberNonZero=0;
+      // and set up mark as char array
+      char * marked = (char *) (index+columnArray->capacity());
+      double * array2 = y->denseVector();
+#ifdef CLP_DEBUG
+      int numberColumns = model->numberColumns();
+      for (i=0;i<numberColumns;i++) {
+        assert(!marked[i]);
+        assert(!array2[i]);
+      }
+#endif
+      for (i=0;i<numberInRowArray;i++) {
+        iRow = whichRow[i]; 
+        double value = pi[i]*scalar;
+        CoinBigIndex j;
+        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array2[iColumn] += value;
+        }
+        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array2[iColumn] -= value;
+        }
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        if (marked[iColumn]) {
+          double value = array2[iColumn];
+          array2[iColumn]=0.0;
+          marked[iColumn]=0;
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+    } else {      
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = mark[i];
+        index[i]=iColumn;
+        array[iColumn]=markVector[iColumn];
+        markVector[iColumn]=0.0;
+      }
+      numberNonZero=numberOriginal;
+      // and set up mark as char array
+      char * marked = (char *) markVector;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        marked[iColumn]=0;
+      }
+      for (i=0;i<numberInRowArray;i++) {
+        iRow = whichRow[i]; 
+        double value = pi[iRow]*scalar;
+        CoinBigIndex j;
+        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array[iColumn] += value;
+        }
+        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
+          int iColumn = column[j];
+          if (!marked[iColumn]) {
+            marked[iColumn]=1;
+            index[numberNonZero++]=iColumn;
+          }
+          array[iColumn] -= value;
+        }
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        marked[iColumn]=0;
+        if (fabs(array[iColumn])>zeroTolerance) {
           index[numberNonZero++]=iColumn;
-        }
-        array[iColumn] += value;
-      }
-      for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
-        int iColumn = column[j];
-        if (!marked[iColumn]) {
-          marked[iColumn]=1;
-          index[numberNonZero++]=iColumn;
-        }
-        array[iColumn] -= value;
-      }
-    }
-    // get rid of tiny values and zero out marked
-    numberOriginal=numberNonZero;
-    numberNonZero=0;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      marked[iColumn]=0;
-      if (fabs(array[iColumn])>zeroTolerance) {
-        index[numberNonZero++]=iColumn;
-      } else {
-        array[iColumn]=0.0;
+        } else {
+          array[iColumn]=0.0;
+        }
       }
     }
   } else if (numberInRowArray==2) {
-    // do by rows when two rows (do longer first)
+    /* do by rows when two rows (do longer first when not packed
+       and shorter first if packed */
     int iRow0 = whichRow[0];
     int iRow1 = whichRow[1];
-    if (startPositive[iRow0+1]-startPositive[iRow0]<
-        startPositive[iRow1+1]-startPositive[iRow1]) {
-      int temp = iRow0;
-      iRow0 = iRow1;
-      iRow1 = temp;
-    }
-    int numberOriginal;
-    int i;
-    numberNonZero=0;
-    double value;
-    value = pi[iRow0]*scalar;
-    CoinBigIndex j;
-    for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
-      int iColumn = column[j];
-      index[numberNonZero++]=iColumn;
-      array[iColumn] = value;
-    }
-    for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
-      int iColumn = column[j];
-      index[numberNonZero++]=iColumn;
-      array[iColumn] = -value;
-    }
-    value = pi[iRow1]*scalar;
-    for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
-      int iColumn = column[j];
-      double value2= array[iColumn];
-      if (value2) {
-        value2 += value;
-      } else {
-        value2 = value;
+    int j;
+    if (packed) {
+      double pi0 = pi[0];
+      double pi1 = pi[1];
+      if (startPositive[iRow0+1]-startPositive[iRow0]>
+          startPositive[iRow1+1]-startPositive[iRow1]) {
+        int temp = iRow0;
+        iRow0 = iRow1;
+        iRow1 = temp;
+        pi0=pi1;
+        pi1=pi[0];
+      }
+      // and set up mark as char array
+      char * marked = (char *) (index+columnArray->capacity());
+      int * lookup = y->getIndices();
+      double value = pi0*scalar;
+      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
+        int iColumn = column[j];
+        array[numberNonZero] = value;
+        marked[iColumn]=1;
+        lookup[iColumn]=numberNonZero;
         index[numberNonZero++]=iColumn;
       }
-      array[iColumn] = value2;
-    }
-    for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
-      int iColumn = column[j];
-      double value2= array[iColumn];
-      if (value2) {
-        value2 -= value;
-      } else {
-        value2 = -value;
+      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
+        int iColumn = column[j];
+        array[numberNonZero] = -value;
+        marked[iColumn]=1;
+        lookup[iColumn]=numberNonZero;
         index[numberNonZero++]=iColumn;
       }
-      array[iColumn] = value2;
-    }
-    // get rid of tiny values and zero out marked
-    numberOriginal=numberNonZero;
-    numberNonZero=0;
-    for (i=0;i<numberOriginal;i++) {
-      int iColumn = index[i];
-      if (fabs(array[iColumn])>zeroTolerance) {
+      int numberOriginal = numberNonZero;
+      value = pi1*scalar;
+      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
+        int iColumn = column[j];
+        if (marked[iColumn]) {
+          int iLookup = lookup[iColumn];
+          array[iLookup] += value;
+        } else {
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero] = value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
+        int iColumn = column[j];
+        if (marked[iColumn]) {
+          int iLookup = lookup[iColumn];
+          array[iLookup] -= value;
+        } else {
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero] = -value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+      // get rid of tiny values and zero out marked
+      int nDelete=0;
+      for (j=0;j<numberOriginal;j++) {
+        int iColumn = index[j];
+        marked[iColumn]=0;
+        if (fabs(array[j])<=zeroTolerance) 
+          nDelete++;
+      }
+      if (nDelete) {
+        numberOriginal=numberNonZero;
+        numberNonZero=0;
+        for (j=0;j<numberOriginal;j++) {
+          int iColumn = index[j];
+          double value = array[j];
+          array[j]=0.0;
+          if (fabs(value)>zeroTolerance) {
+            array[numberNonZero]=value;
+            index[numberNonZero++]=iColumn;
+          }
+        }
+      }
+    } else {
+      if (startPositive[iRow0+1]-startPositive[iRow0]<
+          startPositive[iRow1+1]-startPositive[iRow1]) {
+        int temp = iRow0;
+        iRow0 = iRow1;
+        iRow1 = temp;
+      }
+      int numberOriginal;
+      int i;
+      numberNonZero=0;
+      double value;
+      value = pi[iRow0]*scalar;
+      CoinBigIndex j;
+      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
+        int iColumn = column[j];
         index[numberNonZero++]=iColumn;
-      } else {
-        array[iColumn]=0.0;
+        array[iColumn] = value;
+      }
+      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
+        int iColumn = column[j];
+        index[numberNonZero++]=iColumn;
+        array[iColumn] = -value;
+      }
+      value = pi[iRow1]*scalar;
+      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
+        int iColumn = column[j];
+        double value2= array[iColumn];
+        if (value2) {
+          value2 += value;
+        } else {
+          value2 = value;
+          index[numberNonZero++]=iColumn;
+        }
+        array[iColumn] = value2;
+      }
+      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
+        int iColumn = column[j];
+        double value2= array[iColumn];
+        if (value2) {
+          value2 -= value;
+        } else {
+          value2 = -value;
+          index[numberNonZero++]=iColumn;
+        }
+        array[iColumn] = value2;
+      }
+      // get rid of tiny values and zero out marked
+      numberOriginal=numberNonZero;
+      numberNonZero=0;
+      for (i=0;i<numberOriginal;i++) {
+        int iColumn = index[i];
+        if (fabs(array[iColumn])>zeroTolerance) {
+          index[numberNonZero++]=iColumn;
+        } else {
+          array[iColumn]=0.0;
+        }
       }
     }
@@ -495 +792 @@
     double value;
     iRow = whichRow[0]; 
-    value = pi[iRow]*scalar;
-    if (fabs(value)>zeroTolerance) {
-      CoinBigIndex j;
-      for (j=startPositive[iRow];j<startNegative[iRow];j++) {
-        int iColumn = column[j];
-        index[numberNonZero++]=iColumn;
-        array[iColumn] = value;
-      }
-      for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
-        int iColumn = column[j];
-        index[numberNonZero++]=iColumn;
-        array[iColumn] = -value;
+    CoinBigIndex j;
+    if (packed) {
+      value = pi[0]*scalar;
+      if (fabs(value)>zeroTolerance) {
+        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
+          int iColumn = column[j];
+          array[numberNonZero] = value;
+          index[numberNonZero++]=iColumn;
+        }
+        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
+          int iColumn = column[j];
+          array[numberNonZero] = -value;
+          index[numberNonZero++]=iColumn;
+        }
+      }
+    } else {
+      value = pi[iRow]*scalar;
+      if (fabs(value)>zeroTolerance) {
+        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
+          int iColumn = column[j];
+          array[iColumn] = value;
+          index[numberNonZero++]=iColumn;
+        }
+        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
+          int iColumn = column[j];
+          array[iColumn] = -value;
+          index[numberNonZero++]=iColumn;
+        }
       }
     }
@@ -532 +845 @@
   int numberToDo = y->getNumElements();
   const int * which = y->getIndices();
-  for (jColumn=0;jColumn<numberToDo;jColumn++) {
-    int iColumn = which[jColumn];
-    double value = 0.0;
-    CoinBigIndex j=startPositive_[iColumn];
-    for (;j<startNegative_[iColumn];j++) {
-      int iRow = indices_[j];
-      value += pi[iRow];
-    }
-    for (;j<startPositive_[iColumn+1];j++) {
-      int iRow = indices_[j];
-      value -= pi[iRow];
-    }
-    if (fabs(value)>zeroTolerance) {
-      index[numberNonZero++]=iColumn;
-      array[iColumn]=value;
+  bool packed = rowArray->packedMode();
+  if (packed) {
+    // need to expand pi into y
+    int numberInRowArray = rowArray->getNumElements();
+    assert(y->capacity()>=model->numberRows());
+    double * piOld = pi;
+    pi = y->denseVector();
+    const int * whichRow = rowArray->getIndices();
+    int i;
+    for (i=0;i<numberInRowArray;i++) {
+      int iRow = whichRow[i];
+      pi[iRow]=piOld[i];
+    }
+    // Must line up with y
+    for (jColumn=0;jColumn<numberToDo;jColumn++) {
+      int iColumn = which[jColumn];
+      double value = 0.0;
+      CoinBigIndex j=startPositive_[iColumn];
+      for (;j<startNegative_[iColumn];j++) {
+        int iRow = indices_[j];
+        value += pi[iRow];
+      }
+      for (;j<startPositive_[iColumn+1];j++) {
+        int iRow = indices_[j];
+        value -= pi[iRow];
+      }
+      array[jColumn]=value;
+    }
+    for (i=0;i<numberInRowArray;i++) {
+      int iRow = whichRow[i];
+      pi[iRow]=0.0;
+    }
+  } else {
+    for (jColumn=0;jColumn<numberToDo;jColumn++) {
+      int iColumn = which[jColumn];
+      double value = 0.0;
+      CoinBigIndex j=startPositive_[iColumn];
+      for (;j<startNegative_[iColumn];j++) {
+        int iRow = indices_[j];
+        value += pi[iRow];
+      }
+      for (;j<startPositive_[iColumn+1];j++) {
+        int iRow = indices_[j];
+        value -= pi[iRow];
+      }
+      if (fabs(value)>zeroTolerance) {
+        index[numberNonZero++]=iColumn;
+        array[iColumn]=value;
+      }
     }
   }
@@ -598 +945 @@
   return numberElements;
 }
+/* If element NULL returns number of elements in column part of basis,
+   If not NULL fills in as well */
+CoinBigIndex 
+ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
+                                 const int * whichColumn, 
+                                 int numberBasic,
+                                 int numberColumnBasic,
+                                 int * indexRowU, int * indexColumnU,
+                                 double * elementU) const 
+{
+  int i;
+  CoinBigIndex numberElements=0;
+  if (elementU!=NULL) {
+    assert (columnOrdered_);
+    for (i=0;i<numberColumnBasic;i++) {
+      int iColumn = whichColumn[i];
+      CoinBigIndex j=startPositive_[iColumn];
+      for (;j<startNegative_[iColumn];j++) {
+        int iRow = indices_[j];
+        indexRowU[numberElements]=iRow;
+        indexColumnU[numberElements]=numberBasic;
+        elementU[numberElements++]=1.0;
+      }
+      for (;j<startPositive_[iColumn+1];j++) {
+        int iRow = indices_[j];
+        indexRowU[numberElements]=iRow;
+        indexColumnU[numberElements]=numberBasic;
+        elementU[numberElements++]=-1.0;
+      }
+      numberBasic++;
+    }
+  } else {
+    for (i=0;i<numberColumnBasic;i++) {
+      int iColumn = whichColumn[i];
+      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
+    }
+  }
+  return numberElements;
+}
 /* Unpacks a column into an CoinIndexedvector
-      Note that model is NOT const.  Bounds and objective could
-      be modified if doing column generation */
+ */
 void 
 ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
@@ -615 +1000 @@
     rowArray->add(iRow,-1.0);
   }
+}
+/* Unpacks a column into an CoinIndexedvector
+** in packed foramt
+Note that model is NOT const.  Bounds and objective could
+be modified if doing column generation (just for this variable) */
+void 
+ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
+                            CoinIndexedVector * rowArray,
+                            int iColumn) const
+{
+  int * index = rowArray->getIndices();
+  double * array = rowArray->denseVector();
+  int number = 0;
+  CoinBigIndex j=startPositive_[iColumn];
+  for (;j<startNegative_[iColumn];j++) {
+    int iRow = indices_[j];
+    array[number]=1.0;
+    index[number++]=iRow;
+  }
+  for (;j<startPositive_[iColumn+1];j++) {
+    int iRow = indices_[j];
+    array[number]=-1.0;
+    index[number++]=iRow;
+  }
+  rowArray->setNumElements(number);
+  rowArray->setPackedMode(true);
 }
 /* Adds multiple of a column into an CoinIndexedvector

branches/pre/ClpPrimalColumnDantzig.cpp

@@ -139 +139 @@
 
   for (iSequence=0;iSequence<number;iSequence++) {
-    double value = reducedCost[iSequence];
-    ClpSimplex::Status status = model_->getStatus(iSequence);
-    
-    switch(status) {
+    // check flagged variable
+    if (!model_->flagged(iSequence)) {
+      double value = reducedCost[iSequence];
+      ClpSimplex::Status status = model_->getStatus(iSequence);
       
-    case ClpSimplex::basic:
-    case ClpSimplex::isFixed:
-      break;
-    case ClpSimplex::isFree:
-    case ClpSimplex::superBasic:
-      if (fabs(value)>bestFreeDj) { 
-        bestFreeDj = fabs(value);
-        bestFreeSequence = iSequence;
-      }
-      break;
-    case ClpSimplex::atUpperBound:
-      if (value>bestDj) {
-        bestDj = value;
-        bestSequence = iSequence;
-      }
-      break;
-    case ClpSimplex::atLowerBound:
-      if (value<-bestDj) {
-        bestDj = -value;
-        bestSequence = iSequence;
+      switch(status) {
+
+      case ClpSimplex::basic:
+      case ClpSimplex::isFixed:
+        break;
+      case ClpSimplex::isFree:
+      case ClpSimplex::superBasic:
+        if (fabs(value)>bestFreeDj) { 
+          bestFreeDj = fabs(value);
+          bestFreeSequence = iSequence;
+        }
+        break;
+      case ClpSimplex::atUpperBound:
+        if (value>bestDj) {
+          bestDj = value;
+          bestSequence = iSequence;
+        }
+        break;
+      case ClpSimplex::atLowerBound:
+        if (value<-bestDj) {
+          bestDj = -value;
+          bestSequence = iSequence;
+        }
       }
     }

branches/pre/ClpPrimalColumnSteepest.cpp

@@ -15 +15 @@
 
 // bias for free variables
-
 #define FREE_BIAS 1.0e1
+// Acceptance criteria for free variables
+#define FREE_ACCEPT 1.0e2
 //#############################################################################
 // Constructors / Destructor / Assignment
@@ -227 +228 @@
       }
       if (!model_->nonLinearCost()->lookBothWays()) {
-        
+
         for (j=0;j<number;j++) {
           int iSequence = index[j];
@@ -243 +244 @@
           case ClpSimplex::isFree:
           case ClpSimplex::superBasic:
-            if (fabs(value)>1.0e2*tolerance) {
+            if (fabs(value)>FREE_ACCEPT*tolerance) {
               // we are going to bias towards free (but only if reasonable)
               value *= FREE_BIAS;
@@ -296 +297 @@
           case ClpSimplex::isFree:
           case ClpSimplex::superBasic:
-            if (fabs(value)>1.0e2*tolerance) {
+            if (fabs(value)>FREE_ACCEPT*tolerance) {
               // we are going to bias towards free (but only if reasonable)
               value *= FREE_BIAS;
@@ -377 +378 @@
     case ClpSimplex::isFree:
     case ClpSimplex::superBasic:
-      if (fabs(value)>1.0e2*tolerance) { 
+      if (fabs(value)>FREE_ACCEPT*tolerance) { 
         // we are going to bias towards free (but only if reasonable)
         value *= FREE_BIAS;
@@ -432 +433 @@
         double value = reducedCost[iSequence];
         ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
-        
+
         switch(status) {
           
@@ -481 +482 @@
         double value = reducedCost[iSequence];
         ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
-        
+
         switch(status) {
           
@@ -681 +682 @@
     if (model_->largestDualError()>checkTolerance)
       tolerance *= model_->largestDualError()/checkTolerance;
-  }
+    // But cap
+    tolerance = min(1000.0,tolerance);
+  }
+#ifdef CLP_DEBUG
+  if (model_->numberDualInfeasibilities()==1) 
+    printf("** %g %g %g %x %x %d\n",tolerance,model_->dualTolerance(),
+           model_->largestDualError(),model_,model_->messageHandler(),
+           number);
+#endif
   // stop last one coming immediately
   double saveOutInfeasibility=0.0;
@@ -693 +702 @@
     double value = infeas[iSequence];
     double weight = weights_[iSequence];
-    /*if (model_->numberIterations()%100==0)
-      printf("%d inf %g wt %g\n",iSequence,value,weight);*/
     //weight=1.0;
     if (value>bestDj*weight&&value>tolerance) {
@@ -761 +768 @@
       initializeWeights();
       // create saved weights 
+      delete [] savedWeights_;
       savedWeights_ = new double[numberRows+numberColumns];
       memcpy(savedWeights_,weights_,(numberRows+numberColumns)*
@@ -828 +836 @@
     int number = model_->numberRows() + model_->numberColumns();
     int iSequence;
-
+   
     double * reducedCost = model_->djRegion();
       
@@ -835 +843 @@
         double value = reducedCost[iSequence];
         ClpSimplex::Status status = model_->getStatus(iSequence);
-        
+
         switch(status) {
           
@@ -843 +851 @@
         case ClpSimplex::isFree:
         case ClpSimplex::superBasic:
-          if (fabs(value)>1.0e2*tolerance) { 
+          if (fabs(value)>FREE_ACCEPT*tolerance) { 
             // we are going to bias towards free (but only if reasonable)
             value *= FREE_BIAS;
@@ -867 +875 @@
         double value = reducedCost[iSequence];
         ClpSimplex::Status status = model_->getStatus(iSequence);
-        
+
         switch(status) {
           
@@ -875 +883 @@
         case ClpSimplex::isFree:
         case ClpSimplex::superBasic:
-          if (fabs(value)>1.0e2*tolerance) { 
+          if (fabs(value)>FREE_ACCEPT*tolerance) { 
             // we are going to bias towards free (but only if reasonable)
             value *= FREE_BIAS;
@@ -1158 +1166 @@
   devex_ = 0.0;
 }
+// Returns true if would not find any column
+bool 
+ClpPrimalColumnSteepest::looksOptimal() const
+{
+  //**** THIS MUST MATCH the action coding above
+  double tolerance=model_->currentDualTolerance();
+  // we can't really trust infeasibilities if there is dual error
+  // this coding has to mimic coding in checkDualSolution
+  double error = min(1.0e-3,model_->largestDualError());
+  // allow tolerance at least slightly bigger than standard
+  tolerance = tolerance +  error;
+  if(model_->numberIterations()<model_->lastBadIteration()+200) {
+    // we can't really trust infeasibilities if there is dual error
+    double checkTolerance = 1.0e-8;
+    if (!model_->factorization()->pivots())
+      checkTolerance = 1.0e-6;
+    if (model_->largestDualError()>checkTolerance)
+      tolerance *= model_->largestDualError()/checkTolerance;
+    // But cap
+    tolerance = min(1000.0,tolerance);
+  }
+  int number = model_->numberRows() + model_->numberColumns();
+  int iSequence;
+  
+  double * reducedCost = model_->djRegion();
+  int numberInfeasible=0;
+  if (!model_->nonLinearCost()->lookBothWays()) {
+    for (iSequence=0;iSequence<number;iSequence++) {
+      double value = reducedCost[iSequence];
+      ClpSimplex::Status status = model_->getStatus(iSequence);
+      
+      switch(status) {
+
+      case ClpSimplex::basic:
+      case ClpSimplex::isFixed:
+        break;
+      case ClpSimplex::isFree:
+      case ClpSimplex::superBasic:
+        if (fabs(value)>FREE_ACCEPT*tolerance) 
+          numberInfeasible++;
+        break;
+      case ClpSimplex::atUpperBound:
+        if (value>tolerance) 
+          numberInfeasible++;
+        break;
+      case ClpSimplex::atLowerBound:
+        if (value<-tolerance) 
+          numberInfeasible++;
+      }
+    }
+  } else {
+    ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
+    // can go both ways
+    for (iSequence=0;iSequence<number;iSequence++) {
+      double value = reducedCost[iSequence];
+      ClpSimplex::Status status = model_->getStatus(iSequence);
+      
+      switch(status) {
+
+      case ClpSimplex::basic:
+      case ClpSimplex::isFixed:
+        break;
+      case ClpSimplex::isFree:
+      case ClpSimplex::superBasic:
+        if (fabs(value)>FREE_ACCEPT*tolerance) 
+          numberInfeasible++;
+        break;
+      case ClpSimplex::atUpperBound:
+        if (value>tolerance) {
+          numberInfeasible++;
+        } else {
+          // look other way - change up should be negative
+          value -= nonLinear->changeUpInCost(iSequence);
+          if (value<-tolerance) 
+            numberInfeasible++;
+        }
+        break;
+      case ClpSimplex::atLowerBound:
+        if (value<-tolerance) {
+          numberInfeasible++;
+        } else {
+          // look other way - change down should be positive
+          value -= nonLinear->changeDownInCost(iSequence);
+          if (value>tolerance) 
+            numberInfeasible++;
+        }
+      }
+    }
+  }
+  return numberInfeasible==0;
+}

branches/pre/ClpSimplex.cpp

@@ -104 +104 @@
   numberFake_(0),
   progressFlag_(0),
-  firstFree_(-1)
-
+  firstFree_(-1),
+  incomingInfeasibility_(1.0),
+  allowedInfeasibility_(10.0),
+  progress_(NULL)
 {
   int i;
@@ -125 +127 @@
 }
 
+// Subproblem constructor
+ClpSimplex::ClpSimplex ( const ClpModel * rhs,
+                     int numberRows, const int * whichRow,
+                     int numberColumns, const int * whichColumn,
+                     bool dropNames, bool dropIntegers)
+  : ClpModel(rhs, numberRows, whichRow,
+             numberColumns,whichColumn,dropNames,dropIntegers),
+  columnPrimalInfeasibility_(0.0),
+  rowPrimalInfeasibility_(0.0),
+  columnPrimalSequence_(-2),
+  rowPrimalSequence_(-2), 
+  columnDualInfeasibility_(0.0),
+  rowDualInfeasibility_(0.0),
+  columnDualSequence_(-2),
+  rowDualSequence_(-2),
+  primalToleranceToGetOptimal_(-1.0),
+  remainingDualInfeasibility_(0.0),
+  largeValue_(1.0e15),
+  largestPrimalError_(0.0),
+  largestDualError_(0.0),
+  largestSolutionError_(0.0),
+  dualBound_(1.0e10),
+  alpha_(0.0),
+  theta_(0.0),
+  lowerIn_(0.0),
+  valueIn_(0.0),
+  upperIn_(0.0),
+  dualIn_(0.0),
+  lowerOut_(-1),
+  valueOut_(-1),
+  upperOut_(-1),
+  dualOut_(-1),
+  dualTolerance_(0.0),
+  primalTolerance_(0.0),
+  sumDualInfeasibilities_(0.0),
+  sumPrimalInfeasibilities_(0.0),
+  infeasibilityCost_(1.0e10),
+  sumOfRelaxedDualInfeasibilities_(0.0),
+  sumOfRelaxedPrimalInfeasibilities_(0.0),
+  lower_(NULL),
+  rowLowerWork_(NULL),
+  columnLowerWork_(NULL),
+  upper_(NULL),
+  rowUpperWork_(NULL),
+  columnUpperWork_(NULL),
+  cost_(NULL),
+  rowObjectiveWork_(NULL),
+  objectiveWork_(NULL),
+  sequenceIn_(-1),
+  directionIn_(-1),
+  sequenceOut_(-1),
+  directionOut_(-1),
+  pivotRow_(-1),
+  lastGoodIteration_(-100),
+  dj_(NULL),
+  rowReducedCost_(NULL),
+  reducedCostWork_(NULL),
+  solution_(NULL),
+  rowActivityWork_(NULL),
+  columnActivityWork_(NULL),
+  numberDualInfeasibilities_(0),
+  numberDualInfeasibilitiesWithoutFree_(0),
+  numberPrimalInfeasibilities_(100),
+  numberRefinements_(0),
+  pivotVariable_(NULL),
+  factorization_(NULL),
+  rowScale_(NULL),
+  savedSolution_(NULL),
+  columnScale_(NULL),
+  scalingFlag_(1),
+  numberTimesOptimal_(0),
+  changeMade_(1),
+  algorithm_(0),
+  forceFactorization_(-1),
+  perturbation_(100),
+  nonLinearCost_(NULL),
+  specialOptions_(0),
+  lastBadIteration_(-999999),
+  numberFake_(0),
+  progressFlag_(0),
+  firstFree_(-1),
+  incomingInfeasibility_(1.0),
+  allowedInfeasibility_(10.0),
+  progress_(NULL)
+{
+  int i;
+  for (i=0;i<6;i++) {
+    rowArray_[i]=NULL;
+    columnArray_[i]=NULL;
+  }
+  saveStatus_=NULL;
+  // get an empty factorization so we can set tolerances etc
+  factorization_ = new ClpFactorization();
+  // say Steepest pricing
+  dualRowPivot_ = new ClpDualRowSteepest();
+  // say Steepest pricing
+  primalColumnPivot_ = new ClpPrimalColumnSteepest();
+  solveType_=1; // say simplex based life form
+  
+}
 
 //-----------------------------------------------------------------------------
@@ -140 +242 @@
 }
 int
-ClpSimplex::gutsOfSolution ( const double * rowActivities,
-                             const double * columnActivities,
-                             double * givenDuals,
+ClpSimplex::gutsOfSolution ( double * givenDuals,
                              const double * givenPrimals,
                              bool valuesPass)
@@ -155 +255 @@
     assert(nonLinearCost_);
     int iRow;
-    checkPrimalSolution( rowActivities, columnActivities);
+    checkPrimalSolution( rowActivityWork_, columnActivityWork_);
     // get correct bounds on all variables
-    nonLinearCost_->checkInfeasibilities();
+    nonLinearCost_->checkInfeasibilities(primalTolerance_);
     oldValue = nonLinearCost_->largestInfeasibility();
     save = new double[numberRows_];
@@ -166 +266 @@
   }
   // do work
-  computePrimals(rowActivities, columnActivities);
+  computePrimals(rowActivityWork_, columnActivityWork_);
   // If necessary - override results
   if (givenPrimals) {
@@ -177 +277 @@
     // primal algorithm
     // get correct bounds on all variables
-    nonLinearCost_->checkInfeasibilities();
+    // If  4 bit set - Force outgoing variables to exact bound (primal)
+    if ((specialOptions_&4)==0)
+      nonLinearCost_->checkInfeasibilities(primalTolerance_);
+    else
+      nonLinearCost_->checkInfeasibilities(0.0);
     objectiveModification += nonLinearCost_->changeInCost();
     if (nonLinearCost_->numberInfeasibilities())
@@ -191 +295 @@
 #endif
     int numberOut=0;
-    if (oldValue<1.0
-        &&nonLinearCost_->largestInfeasibility()>1.1*oldValue+1.0e-4||
+    if (oldValue<incomingInfeasibility_
+        &&nonLinearCost_->largestInfeasibility()>
+        max(incomingInfeasibility_,allowedInfeasibility_)||
         largestPrimalError_>1.0e-3) {
+      printf("Original largest infeas %g, now %g, primalError %g\n",
+             oldValue,nonLinearCost_->largestInfeasibility(),
+             largestPrimalError_);
       // throw out up to 1000 structurals
       int iRow;
@@ -200 +308 @@
       for (iRow=0;iRow<numberRows_;iRow++) {
         int iPivot=pivotVariable_[iRow];
-        double difference = fabs(solution_[iPivot]=save[iRow]);
+        double difference = fabs(solution_[iPivot]-save[iRow]);
         solution_[iPivot]=save[iRow];
         save[iRow]=difference;
@@ -234 +342 @@
 
   // now check solutions
-  checkPrimalSolution( rowActivities, columnActivities);
+  checkPrimalSolution( rowActivityWork_, columnActivityWork_);
   objectiveValue_ += objectiveModification;
   checkDualSolution();
@@ -243 +351 @@
       <<largestDualError_
       <<CoinMessageEol;
+  // Switch off false values pass indicator
+  if (!valuesPass&&algorithm_>0)
+    firstFree_ = -1;
   return 0;
 }
@@ -253 +364 @@
   CoinIndexedVector  * workSpace = rowArray_[0];
 
-  double * array = new double [numberRows_+1]; // +1 for network
-  double * save = new double [numberRows_+1];
-  double * previous = new double [numberRows_+1];
+  CoinIndexedVector arrayVector;
+  arrayVector.reserve(numberRows_+1);
+  CoinIndexedVector previousVector;
+  previousVector.reserve(numberRows_+1);
 
   // accumulate non basic stuff 
-  ClpFillN(array,numberRows_,0.0);
 
   int iRow;
@@ -272 +383 @@
     solution_[iPivot] = 0.0;
   }
+  double * array = arrayVector.denseVector();
   times(-1.0,columnActivityWork_,array);
 
+  int * index = arrayVector.getIndices();
+  int number=0;
   for (iRow=0;iRow<numberRows_;iRow++) {
-    array[iRow] += rowActivityWork_[iRow];
-  }
+    double value = array[iRow] + rowActivityWork_[iRow];
+    if (value) {
+      array[iRow]=value;
+      index[number++]=iRow;
+    } else {
+      array[iRow]=0.0;
+    }
+  }
+
+  arrayVector.setNumElements(number);
 
   // Ftran adjusted RHS and iterate to improve accuracy
@@ -282 +404 @@
   int iRefine;
   double * work = workSpace->denseVector();
-  factorization_->updateColumn(workSpace,array);
-
+  CoinIndexedVector * thisVector = &arrayVector;
+  CoinIndexedVector * lastVector = &previousVector;
+  factorization_->updateColumn(workSpace,thisVector);
+  bool goodSolution=true;
   for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
-    // save in case we want to restore
-    ClpDisjointCopyN ( array, numberRows_ , save);
-
+
+    int numberIn = thisVector->getNumElements();
+    int * indexIn = thisVector->getIndices();
+    double * arrayIn = thisVector->denseVector();
+    // put solution in correct place
+    int j;
+    for (j=0;j<numberIn;j++) {
+      iRow = indexIn[j];
+      int iPivot=pivotVariable_[iRow];
+      solution_[iPivot] = arrayIn[iRow];
+    }
+    // check Ax == b  (for all)
+    times(-1.0,columnActivityWork_,work);
+    largestPrimalError_=0.0;
+    double multiplier = 131072.0;
+    for (iRow=0;iRow<numberRows_;iRow++) {
+      double value = work[iRow] + rowActivityWork_[iRow];
+      work[iRow] = value*multiplier;
+      if (fabs(value)>largestPrimalError_) {
+        largestPrimalError_=fabs(value);
+      }
+    }
+    if (largestPrimalError_>=lastError) {
+      // restore
+      CoinIndexedVector * temp = thisVector;
+      thisVector = lastVector;
+      lastVector=temp;
+      goodSolution=false;
+      break;
+    }
+    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-10) {
+      // try and make better
+      // save this
+      CoinIndexedVector * temp = thisVector;
+      thisVector = lastVector;
+      lastVector=temp;
+      int * indexOut = thisVector->getIndices();
+      int number=0;
+      array = thisVector->denseVector();
+      thisVector->clear();
+      for (iRow=0;iRow<numberRows_;iRow++) {
+        double value = work[iRow];
+        if (value) {
+          array[iRow]=value;
+          indexOut[number++]=iRow;
+          work[iRow]=0.0;
+        }
+      }
+      thisVector->setNumElements(number);
+      lastError=largestPrimalError_;
+      factorization_->updateColumn(workSpace,thisVector);
+      multiplier = 1.0/multiplier;
+      double * previous = lastVector->denseVector();
+      number=0;
+      for (iRow=0;iRow<numberRows_;iRow++) {
+        double value = previous[iRow] + multiplier*array[iRow];
+        if (value) {
+          array[iRow]=value;
+          indexOut[number++]=iRow;
+        } else {
+          array[iRow]=0.0;
+        }
+      }
+      thisVector->setNumElements(number);
+    } else {
+      break;
+    }
+  }
+
+  // solution as accurate as we are going to get
+  ClpFillN(work,numberRows_,0.0);
+  if (!goodSolution) {
+    array = thisVector->denseVector();
     // put solution in correct place
     for (iRow=0;iRow<numberRows_;iRow++) {
@@ -293 +487 @@
       solution_[iPivot] = array[iRow];
     }
-    // check Ax == b  (for all)
-    times(-1.0,columnActivityWork_,work);
-    for (iRow=0;iRow<numberRows_;iRow++) {
-      work[iRow] += rowActivityWork_[iRow];
-    }
-
-    largestPrimalError_=0.0;
-    for (iRow=0;iRow<numberRows_;iRow++) {
-      if (fabs(work[iRow])>largestPrimalError_) {
-        largestPrimalError_=fabs(work[iRow]);
-      }
-      //work[iRow] -= save[iRow];
-    }
-    if (largestPrimalError_>=lastError) {
-      // restore
-      double * temp = array;
-      array = previous;
-      previous=temp;
-      break;
-    }
-    if (iRefine<numberRefinements_&&largestPrimalError_>1.0e-12) {
-      // try and make better
-      // save this
-      double * temp = array;
-      array = previous;
-      previous=temp;
-      double multiplier = 131072.0;
-      for (iRow=0;iRow<numberRows_;iRow++) {
-        array[iRow] = multiplier*work[iRow];
-        work[iRow]=0.0;
-      }
-      lastError=largestPrimalError_;
-      factorization_->updateColumn(workSpace,array);
-      multiplier = 1.0/multiplier;
-      for (iRow=0;iRow<numberRows_;iRow++) {
-        array[iRow] = previous[iRow] + multiplier*array[iRow];
-        work[iRow]=0.0;
-      }
-    }
-  }
-
-  // solution as accurate as we are going to get
-  ClpFillN(work,numberRows_,0.0);
-  // put solution in correct place
-  for (iRow=0;iRow<numberRows_;iRow++) {
-    int iPivot=pivotVariable_[iRow];
-    solution_[iPivot] = array[iRow];
-  }
-  delete [] previous;
-  delete [] array;
-  delete [] save;
+  }
 }
 // now dual side
@@ -352 +496 @@
   CoinIndexedVector  * workSpace = rowArray_[0];
 
-  double * array = new double [numberRows_+1]; // +1 for network
-  double * save = new double [numberRows_+1];
-  double * previous = new double [numberRows_+1];
+  CoinIndexedVector arrayVector;
+  arrayVector.reserve(numberRows_+1);
+  CoinIndexedVector previousVector;
+  previousVector.reserve(numberRows_+1);
+
 
   int iRow;
@@ -360 +506 @@
   workSpace->checkClear();
 #endif
+  double * array = arrayVector.denseVector();
+  int * index = arrayVector.getIndices();
+  int number=0;
   if (!givenDjs) {
     for (iRow=0;iRow<numberRows_;iRow++) {
       int iPivot=pivotVariable_[iRow];
-      array[iRow]=cost_[iPivot];
+      double value = cost_[iPivot];
+      if (value) {
+        array[iRow]=value;
+        index[number++]=iRow;
+      }
     }
   } else {
@@ -372 +525 @@
       if (!pivoted(iPivot))
         givenDjs[iPivot]=0.0;
-      array[iRow]=cost_[iPivot]-givenDjs[iPivot];
-    }
-  }
-  ClpDisjointCopyN ( array, numberRows_ , save);
+      double value =cost_[iPivot]-givenDjs[iPivot];
+      if (value) {
+        array[iRow]=value;
+        index[number++]=iRow;
+      }
+    }
+  }
+  arrayVector.setNumElements(number);
 
   // Btran basic costs and get as accurate as possible
@@ -381 +538 @@
   int iRefine;
   double * work = workSpace->denseVector();
-  factorization_->updateColumnTranspose(workSpace,array);
+  CoinIndexedVector * thisVector = &arrayVector;
+  CoinIndexedVector * lastVector = &previousVector;
+  factorization_->updateColumnTranspose(workSpace,thisVector);
+
   for (iRefine=0;iRefine<numberRefinements_+1;iRefine++) {
     // check basic reduced costs zero
@@ -391 +551 @@
       for (iRow=0;iRow<numberRows_;iRow++) {
         int iPivot=pivotVariable_[iRow];
+        double value;
         if (iPivot>=numberColumns_) {
           // slack
-          work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
+          value = rowObjectiveWork_[iPivot-numberColumns_]
             + array[iPivot-numberColumns_];
         } else {
           // column
-        work[iRow] = reducedCostWork_[iPivot];
-        }
-        if (fabs(work[iRow])>largestDualError_) {
-          largestDualError_=fabs(work[iRow]);
+          value = reducedCostWork_[iPivot];
+        }
+        work[iRow]=value;
+        if (fabs(value)>largestDualError_) {
+          largestDualError_=fabs(value);
         }
       }
@@ -424 +586 @@
     if (largestDualError_>=lastError) {
       // restore
-      double * temp = array;
-      array = previous;
-      previous=temp;
+      CoinIndexedVector * temp = thisVector;
+      thisVector = lastVector;
+      lastVector=temp;
       break;
     }
-    if (iRefine<numberRefinements_&&largestDualError_>1.0e-20
+    if (iRefine<numberRefinements_&&largestDualError_>1.0e-10
         &&!givenDjs) {
       // try and make better
       // save this
-      double * temp = array;
-      array = previous;
-      previous=temp;
+      CoinIndexedVector * temp = thisVector;
+      thisVector = lastVector;
+      lastVector=temp;
+      int * indexOut = thisVector->getIndices();
+      int number=0;
+      array = thisVector->denseVector();
+      thisVector->clear();
       double multiplier = 131072.0;
       for (iRow=0;iRow<numberRows_;iRow++) {
-        array[iRow] = multiplier*work[iRow];
+        double value = multiplier*work[iRow];
+        if (value) {
+          array[iRow]=value;
+          indexOut[number++]=iRow;
+          work[iRow]=0.0;
+        }
         work[iRow]=0.0;
       }
+      thisVector->setNumElements(number);
       lastError=largestDualError_;
-      factorization_->updateColumnTranspose(workSpace,array);
+      factorization_->updateColumnTranspose(workSpace,thisVector);
       multiplier = 1.0/multiplier;
+      double * previous = lastVector->denseVector();
+      number=0;
       for (iRow=0;iRow<numberRows_;iRow++) {
-        array[iRow] = previous[iRow] + multiplier*array[iRow];
-        work[iRow]=0.0;
-      }
+        double value = previous[iRow] + multiplier*array[iRow];
+        if (value) {
+          array[iRow]=value;
+          indexOut[number++]=iRow;
+        } else {
+          array[iRow]=0.0;
+        }
+      }
+      thisVector->setNumElements(number);
     } else {
       break;
@@ -454 +634 @@
   ClpFillN(work,numberRows_,0.0);
   // now look at dual solution
+  array = thisVector->denseVector();
   for (iRow=0;iRow<numberRows_;iRow++) {
     // slack
@@ -469 +650 @@
   }
 
-  delete [] previous;
-  delete [] array;
-  delete [] save;
 }
 /* Given an existing factorization computes and checks 
@@ -483 +661 @@
     createRim(7+8+16+32);
     // do work
-    gutsOfSolution ( rowActivities, columnActivities,NULL,NULL);
+    gutsOfSolution ( NULL,NULL);
     // release extra memory
     deleteRim(0);
@@ -519 +697 @@
 /* Factorizes using current basis.  
       solveType - 1 iterating, 0 initial, -1 external 
+      - 2 then iterating but can throw out of basis
       If 10 added then in primal values pass
 */
@@ -528 +707 @@
   int iRow,iColumn;
   int totalSlacks=numberRows_;
+  if (!status_)
+    createStatus();
 
   bool valuesPass=false;
@@ -535 +716 @@
   }
 #ifdef CLP_DEBUG
-  if (solveType==1) {
+  if (solveType>0) {
     int numberFreeIn=0,numberFreeOut=0;
     double biggestDj=0.0;
     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
       switch(getColumnStatus(iColumn)) {
-        
+
       case basic:
         if (columnLower_[iColumn]<-largeValue_
@@ -560 +741 @@
   }
 #endif
-  if (!solveType) {
+  if (solveType<=0) {
+    // Make sure everything is clean
+    for (iRow=0;iRow<numberRows_;iRow++) {
+      if(getRowStatus(iRow)==isFixed) {
+        // double check fixed
+        if (rowUpperWork_[iRow]>rowLowerWork_[iRow])
+          setRowStatus(iRow,atLowerBound);
+      } else if (getRowStatus(iRow)==isFree) {
+        // may not be free after all
+        if (rowLowerWork_[iRow]>-largeValue_||rowUpperWork_[iRow]<largeValue_)
+          setRowStatus(iRow,superBasic);
+      }
+    }
+    for (iColumn=0;iColumn<numberColumns_;iColumn++) {
+      if(getColumnStatus(iColumn)==isFixed) {
+        // double check fixed
+        if (columnUpperWork_[iColumn]>columnLowerWork_[iColumn])
+          setColumnStatus(iColumn,atLowerBound);
+      } else if (getColumnStatus(iColumn)==isFree) {
+        // may not be free after all
+        if (columnLowerWork_[iColumn]>-largeValue_||columnUpperWork_[iColumn]<largeValue_)
+          setColumnStatus(iColumn,superBasic);
+      }
+    }
     if (!valuesPass) {
       // not values pass so set to bounds
@@ -609 +813 @@
             break;
           case isFree:
-            // may not be free after all
-            if (rowLowerWork_[iRow]<-largeValue_&&rowUpperWork_[iRow]>largeValue_)
               break;
             // not really free - fall through to superbasic
@@ -682 +884 @@
             }
             break;
+          case isFixed:
           case atLowerBound:
-          case ClpSimplex::isFixed:
             columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
             if (columnActivityWork_[iColumn]<-largeValue_) {
@@ -696 +898 @@
             break;
           case isFree:
-            // may not be free after all
-            if (columnLowerWork_[iColumn]<-largeValue_&&columnUpperWork_[iColumn]>largeValue_)
               break;
             // not really free - fall through to superbasic
@@ -881 +1081 @@
     }
   }
+#if 0
   int status=-99;
   int * rowIsBasic = new int[numberRows_];
@@ -907 +1108 @@
     assert (numberBasic<=numberRows_);
     while (status==-99) {
-      status =  factorization_->factorize(this,matrix_,
+      status =  factorization_->factorizeOld(this,matrix_,
                                           numberRows_,numberColumns_,
                                           rowIsBasic, columnIsBasic,
@@ -917 +1118 @@
     }
     if (!status) {
-  /// See whether to redo pivot order
-      if (factorization_->needToReorder()||!numberIterations_) {
+      // See whether to redo pivot order
+      if (factorization_->needToReorder()||
+          progress_->lastIterationNumber(0)<0) {
         // do pivot information
         for (i=0;i<numberRows_;i++) {
@@ -1018 +1220 @@
   delete [] rowIsBasic;
   delete [] columnIsBasic;
+#else
+  int status = factorization_->factorize(this, solveType,valuesPass);
+#endif
   if (status) {
     handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)
@@ -1108 +1313 @@
   if (hitMaximumIterations())
     return 2;
+#if 0
+  // check for small cycles
+  int cycle=progress_->cycle(sequenceIn_,sequenceOut_,
+                            directionIn_,directionOut_);
+  if (cycle>0) {
+    printf("Cycle of %d\n",cycle);
+    if (factorization_->pivots()>cycle)
+      forceFactorization_=cycle-1;
+    else
+      forceFactorization_=1;
+    return 1;
+  }
+#endif
   // only time to re-factorize if one before real time
   // this is so user won't be surprised that maximumPivots has exact meaning
@@ -1116 +1334 @@
         factorization_->pivots()==forceFactorization_) {
       // relax
-      forceFactorization_ = (1+3*forceFactorization_)/2;
+      forceFactorization_ = (3+5*forceFactorization_)/4;
       if (forceFactorization_>factorization_->maximumPivots())
         forceFactorization_ = -1; //off
@@ -1202 +1420 @@
   numberFake_(0),
   progressFlag_(0),
-  firstFree_(-1)
+  firstFree_(-1),
+  incomingInfeasibility_(1.0),
+  allowedInfeasibility_(10.0),
+  progress_(NULL)
 {
   int i;
@@ -1296 +1517 @@
   numberFake_(0),
   progressFlag_(0),
-  firstFree_(-1)
+  firstFree_(-1),
+  incomingInfeasibility_(1.0),
+  allowedInfeasibility_(10.0),
+  progress_(NULL)
 {
   int i;
@@ -1308 +1532 @@
   // say Steepest pricing
   dualRowPivot_ = new ClpDualRowSteepest();
-  // say Dantzig pricing
-  primalColumnPivot_ = new ClpPrimalColumnDantzig();
+  // say Steepest pricing
+  primalColumnPivot_ = new ClpPrimalColumnSteepest();
   solveType_=1; // say simplex based life form
   
@@ -1336 +1560 @@
   rowUpperWork_ = upper_+numberColumns_;
   columnUpperWork_ = upper_;
+  //cost_ = ClpCopyOfArray(rhs.cost_,2*(numberColumns_+numberRows_));
   cost_ = ClpCopyOfArray(rhs.cost_,(numberColumns_+numberRows_));
   objectiveWork_ = cost_;
@@ -1429 +1654 @@
   progressFlag_ = rhs.progressFlag_;
   firstFree_ = rhs.firstFree_;
+  incomingInfeasibility_ = rhs.incomingInfeasibility_;
+  allowedInfeasibility_ = rhs.allowedInfeasibility_;
+  if (rhs.progress_)
+    progress_ = new ClpSimplexProgress(*rhs.progress_);
+  else
+    progress_=NULL;
   sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
   sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
@@ -1492 +1723 @@
     delete primalColumnPivot_;
     primalColumnPivot_ = NULL;
+    delete progress_;
+    progress_=NULL;
   } else {
     // delete any size information in methods
@@ -1512 +1745 @@
 
   objectiveValue_ = 0.0;
-  firstFree_ = -1;
   // now look at primal solution
   columnPrimalInfeasibility_=0.0;
@@ -1590 +1822 @@
   rowDualInfeasibility_=0.0;
   rowDualSequence_=-1;
-  int firstFree = -1;
+  int firstFreePrimal = -1;
+  int firstFreeDual = -1;
+  int numberSuperBasicWithDj=0;
   primalToleranceToGetOptimal_=max(rowPrimalInfeasibility_,
                                    columnPrimalInfeasibility_);
@@ -1602 +1836 @@
 
   for (iColumn=0;iColumn<numberColumns_;iColumn++) {
-    if (getColumnStatus(iColumn) != basic) {
+    if (getColumnStatus(iColumn) != basic&&!flagged(iColumn)) {
       // not basic
       double distanceUp = columnUpperWork_[iColumn]-
@@ -1610 +1844 @@
       if (distanceUp>primalTolerance_) {
         double value = reducedCostWork_[iColumn];
-        if (!flagged(iColumn)) {
-          // Check if "free"
-          if (firstFree<0&&distanceDown>primalTolerance_) {
-            if (algorithm_>0) {
-              // primal
-              firstFree = iColumn;
-            } else if (fabs(value)>1.0e2*relaxedTolerance) {
-              // dual with dj
-              firstFree = iColumn;
-            }
+        // Check if "free"
+        if (distanceDown>primalTolerance_) {
+          if (fabs(value)>1.0e2*relaxedTolerance) {
+            numberSuperBasicWithDj++;
+            if (firstFreeDual<0)
+              firstFreeDual = iColumn;
           }
+          if (firstFreePrimal<0)
+            firstFreePrimal = iColumn;
         }
         // should not be negative
@@ -1683 +1915 @@
   }
   for (iRow=0;iRow<numberRows_;iRow++) {
-    if (getRowStatus(iRow) != basic) {
+    if (getRowStatus(iRow) != basic&&!flagged(iRow+numberColumns_)) {
       // not basic
       double distanceUp = rowUpperWork_[iRow]-rowActivityWork_[iRow];
@@ -1689 +1921 @@
       if (distanceUp>primalTolerance_) {
         double value = rowReducedCost_[iRow];
-        if (!flagged(iRow+numberColumns_)) {
-          // Check if "free"
-          if (firstFree<0&&distanceDown>primalTolerance_) {
-            if (algorithm_>0) {
-              // primal
-              firstFree = iRow+numberColumns_;
-            } else if (fabs(value)>1.0e2*relaxedTolerance) {
-              // dual with dj
-              firstFree = iRow+numberColumns_;
-            }
+        // Check if "free"
+        if (distanceDown>primalTolerance_) {
+          if (fabs(value)>1.0e2*relaxedTolerance) {
+            numberSuperBasicWithDj++;
+            if (firstFreeDual<0)
+              firstFreeDual = iRow+numberColumns_;
           }
+          if (firstFreePrimal<0)
+            firstFreePrimal = iRow+numberColumns_;
         }
         // should not be negative
@@ -1751 +1981 @@
     }
   }
-  if (firstFree>=0) {
-    if (algorithm_<0||!numberIterations_)
-      firstFree_=firstFree;
+  if (algorithm_<0&&firstFreeDual>=0) {
+    // dual
+    firstFree_ = firstFreeDual;
+  } else if (numberSuperBasicWithDj||
+             (progress_&&progress_->lastIterationNumber(0)<=0)) {
+    firstFree_=firstFreePrimal;
   }
 }
@@ -1783 +2016 @@
   }
 }
+/*
+  Unpacks one column of the matrix into indexed array 
+*/
+void 
+ClpSimplex::unpackPacked(CoinIndexedVector * rowArray) 
+{
+  rowArray->clear();
+  if (sequenceIn_>=numberColumns_) {
+    //slack
+    int * index = rowArray->getIndices();
+    double * array = rowArray->denseVector();
+    array[0]=-1.0;
+    index[0]=sequenceIn_-numberColumns_;
+    rowArray->setNumElements(1);
+    rowArray->setPackedMode(true);
+  } else {
+    // column
+    matrix_->unpackPacked(this,rowArray,sequenceIn_);
+  }
+}
+void 
+ClpSimplex::unpackPacked(CoinIndexedVector * rowArray,int sequence)
+{
+  rowArray->clear();
+  if (sequence>=numberColumns_) {
+    //slack
+    int * index = rowArray->getIndices();
+    double * array = rowArray->denseVector();
+    array[0]=-1.0;
+    index[0]=sequence-numberColumns_;
+    rowArray->setNumElements(1);
+    rowArray->setPackedMode(true);
+  } else {
+    // column
+    matrix_->unpackPacked(this,rowArray,sequence);
+  }
+}
 bool
 ClpSimplex::createRim(int what,bool makeRowCopy)
 {
   bool goodMatrix=true;
+  int saveLevel=handler_->logLevel();
+  if (problemStatus_==10)
+    handler_->setLogLevel(0); // switch off messages
   int i;
   if ((what&(16+32))!=0) {
@@ -1845 +2118 @@
   if ((what&4)!=0) {
     delete [] cost_;
-    cost_ = new double[numberColumns_+numberRows_];
+    // extra copy with original costs
+    int nTotal = numberRows_+numberColumns_;
+    //cost_ = new double[2*nTotal];
+    cost_ = new double[nTotal];
     objectiveWork_ = cost_;
     rowObjectiveWork_ = cost_+numberColumns_;
@@ -1853 +2129 @@
     else
       memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
+    // and initialize changes to zero
+    //memset(cost_+nTotal,0,nTotal*sizeof(double));
   }
   if ((what&1)!=0) {
@@ -1915 +2193 @@
       if (columnUpperWork_[i]<1.0e50)
         columnUpperWork_[i] *= multiplier;
+      
     }
     for (i=0;i<numberRows_;i++) {
@@ -1971 +2250 @@
     }    
   }
+  double primalTolerance=dblParam_[ClpPrimalTolerance];
+  if ((what&1)!=0) {
+    // fix any variables with tiny gaps
+    for (i=0;i<numberColumns_;i++) {
+      if (columnUpperWork_[i]-columnLowerWork_[i]<=primalTolerance) {
+        if (columnLowerWork_[i]>=0.0) {
+          columnUpperWork_[i] = columnLowerWork_[i];
+        } else if (columnUpperWork_[i]<=0.0) {
+          columnLowerWork_[i] = columnUpperWork_[i];
+        } else {
+          columnUpperWork_[i] = 0.0;
+          columnLowerWork_[i] = 0.0;
+        }
+      }
+    }
+    for (i=0;i<numberRows_;i++) {
+      if (rowUpperWork_[i]-rowLowerWork_[i]<=primalTolerance) {
+        if (rowLowerWork_[i]>=0.0) {
+          rowUpperWork_[i] = rowLowerWork_[i];
+        } else if (rowUpperWork_[i]<=0.0) {
+          rowLowerWork_[i] = rowUpperWork_[i];
+        } else {
+          rowUpperWork_[i] = 0.0;
+          rowLowerWork_[i] = 0.0;
+        }
+      }
+    }
+  }
+  if (problemStatus_==10) {
+    problemStatus_=-1;
+    handler_->setLogLevel(saveLevel); // switch back messages
+  }
   return goodMatrix;
 }
@@ -2037 +2348 @@
     infeasibilityCost_=value;
 }
-void ClpSimplex::setnumberRefinements( int value) 
+void ClpSimplex::setNumberRefinements( int value) 
 {
   if (value>=0&&value<10)
@@ -2179 +2490 @@
     }
   }
-#define MAXPASS 10  
+#define MAXPASS 10
 
   // Loop round seeing if we can tighten bounds
   // Would be faster to have a stack of possible rows
   // and we put altered rows back on stack
-
-  while(numberChanged) {
+  int numberCheck=-1;
+  while(numberChanged>numberCheck) {
 
     numberChanged = 0; // Bounds tightened this pass
@@ -2205 +2516 @@
         CoinBigIndex rEnd = rowStart[iRow]+rowLength[iRow];
         CoinBigIndex j;
-
         // Compute possible lower and upper ranges
-
+      
         for (j = rStart; j < rEnd; ++j) {
           double value=element[j];
@@ -2213 +2523 @@
           if (value > 0.0) {
             if (columnUpper_[iColumn] >= large) {
-              maximumUp = 1e31;
               ++infiniteUpper;
             } else {
@@ -2219 +2528 @@
             }
             if (columnLower_[iColumn] <= -large) {
-              maximumDown = -1e31;
               ++infiniteLower;
             } else {
@@ -2226 +2534 @@
           } else if (value<0.0) {
             if (columnUpper_[iColumn] >= large) {
-              maximumDown = -1e31;
               ++infiniteLower;
             } else {
@@ -2232 +2539 @@
             }
             if (columnLower_[iColumn] <= -large) {
-              maximumUp = 1e31;
               ++infiniteUpper;
             } else {
@@ -2239 +2545 @@
           }
         }
-        if (maximumUp <= rowUpper_[iRow] + tolerance && 
-            maximumDown >= rowLower_[iRow] - tolerance) {
-
+        double maxUp = maximumUp+infiniteUpper*1.0e31;
+        double maxDown = maximumDown-infiniteLower*1.0e31;
+        if (maxUp <= rowUpper_[iRow] + tolerance && 
+            maxDown >= rowLower_[iRow] - tolerance) {
+          
           // Row is redundant - make totally free
-          rowLower_[iRow]=-COIN_DBL_MAX;
-          rowUpper_[iRow]=COIN_DBL_MAX;
+          // NO - can't do this for postsolve
+          // rowLower_[iRow]=-COIN_DBL_MAX;
+          // rowUpper_[iRow]=COIN_DBL_MAX;
+          //printf("Redundant row in presolveX %d\n",iRow);
 
         } else {
-          if (maximumUp < rowLower_[iRow] -tolerance ||
-              maximumDown > rowUpper_[iRow]+tolerance) {
+          if (maxUp < rowLower_[iRow] -tolerance ||
+              maxDown > rowUpper_[iRow]+tolerance) {
             // problem is infeasible - exit at once
             numberInfeasible++;
-            if (handler_->logLevel()>3)
-              printf("a %g %g %g %g\n",maximumUp,rowLower_[iRow],
-                     maximumDown,rowUpper_[iRow]);
             break;
           }
-
-          if (infiniteUpper == 0 && rowLower_[iRow] > -large) {
-            for (j = rStart; j < rEnd; ++j) {
-              double value=element[j];
-              iColumn = column[j];
-              if (value > 0.0) {
-                if (columnUpper_[iColumn] < large) {
-                  newBound = columnUpper_[iColumn] + 
-                    (rowLower_[iRow] - maximumUp) / value;
-                  if (newBound > columnLower_[iColumn] + 1.0e-12) {
-                    // Tighten the lower bound 
-
-                    columnLower_[iColumn] = newBound;
-                    ++numberChanged;
-
-                    // check infeasible (relaxed)
-                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
-                        -100.0*tolerance) {
-                      numberInfeasible++;
-                      if (handler_->logLevel()>3)
-                        printf("b %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
-                    }
-                    infiniteLower=1; // skip looking at other bound
+          double lower = rowLower_[iRow];
+          double upper = rowUpper_[iRow];
+          for (j = rStart; j < rEnd; ++j) {
+            double value=element[j];
+            iColumn = column[j];
+            double nowLower = columnLower_[iColumn];
+            double nowUpper = columnUpper_[iColumn];
+            if (value > 0.0) {
+              // positive value
+              if (lower>-large) {
+                if (!infiniteUpper) {
+                  assert(nowUpper < large);
+                  newBound = nowUpper + 
+                    (lower - maximumUp) / value;
+                } else if (infiniteUpper==1&&nowUpper>large) {
+                  newBound = (lower -maximumUp) / value;
+                } else {
+                  newBound = -COIN_DBL_MAX;
+                }
+                if (newBound > nowLower + 1.0e-12) {
+                  // Tighten the lower bound 
+                  columnLower_[iColumn] = newBound;
+                  numberChanged++;
+                  // check infeasible (relaxed)
+                  if (nowUpper - newBound < 
+                      -100.0*tolerance) {
+                    numberInfeasible++;
                   }
+                  // adjust
+                  double now;
+                  if (nowLower<-large) {
+                    now=0.0;
+                    infiniteLower--;
+                  } else {
+                    now = nowLower;
+                  }
+                  maximumDown += (newBound-now) * value;
+                  nowLower = newBound;
                 }
-              } else {
-                if (columnLower_[iColumn] > -large) {
-                  newBound = columnLower_[iColumn] + 
-                    (rowLower_[iRow] - maximumUp) / value;
-                  if (newBound < columnUpper_[iColumn] - 1.0e-12) {
-                    // Tighten the upper bound 
-
-                    columnUpper_[iColumn] = newBound;
-                    ++numberChanged;
-
-                    // check infeasible (relaxed)
-                    if (columnUpper_[iColumn] - columnLower_[iColumn] < 
-                        -100.0*tolerance) {
-                      numberInfeasible++;
-                      if (handler_->logLevel()>3)
-                        printf("c %g %g \n",columnLower_[iColumn],columnUpper_[iColumn]);
-                    }
-                    infiniteLower=1; // skip looking at other bound
+              } 
+              if (upper <large) {
+                if (!infiniteLower) {
+                  assert(nowLower >- large);
+                  newBound = nowLower + 
+                    (upper - maximumDown) / value;
+                } else if (infiniteLower==1&&nowLower<-large) {
+                  newBound =   (upper - maximumDown) / value;
+                } else {
+                  newBound = COIN_DBL_MAX;
+                }
+                if (newBound < nowUpper - 1.0e-12) {
+                  // Tighten the upper bound 
+                  columnUpper_[iColumn] = newBound;
+                  numberChanged++;
+                  // check infeasible (relaxed)
+                  if (newBound - nowLower < 
+                      -100.0*tolerance) {
+                    numberInfeasible++;
                   }
+                  // adjust 
+                  double now;
+                  if (nowUpper>large) {
+                    now=0.0;
+                    infiniteUpper--;
+                  } else {
+                    now = nowUpper;
+                  }
+                  maximumUp += (newBound-now) * value;
+                  nowUpper = newBound;
+                }
+              }
+            } else {
+              // negative value
+              if (lower>-large) {
+                if (!infiniteUpper) {
+                  assert(nowLower >- large);
+                  newBound = nowLower + 
+                    (lower - maximumUp) / value;
+                } else if (infiniteUpper==1&&nowLower<-large) {
+                  newBound = (lower -maximumUp) / value;
+                } else {
+                  newBound = COIN_DBL_MAX;
+                }
+                if (newBound < nowUpper - 1.0e-12) {
+                  // Tighten the upper bound 
+                  columnUpper_[iColumn] = newBound;
+                  numberChanged++;
+                  // check infeasible (relaxed)
+                  if (newBound - nowLower < 
+                      -100.0*tolerance) {
+                    numberInfeasible++;
+                  }
+                  // adjust
+                  double now;
+                  if (nowUpper>large) {
+                    now=0.0;
+                    infiniteLower--;
+                  } else {
+                    now = nowUpper;
+                  }
+                  maximumDown += (newBound-now) * value;
+                  nowUpper = newBound;
+                }
+              }
+              if (upper <large) {
+                if (!infiniteLower) {
+                  assert(nowUpper < large);
+                  newBound = nowUpper + 
+                    (upper - maximumDown) / value;
+                } else if (infiniteLower==1&&nowUpper>large) {
+                  newBound =   (upper - maximumDown) / value;
+                } else {
+                  newBound = -COIN_DBL_MAX;
+                }
+                if (newBound > nowLower + 1.0e-12) {
+                  // Tighten the lower bound 
+                  columnLower_[iColumn] = newBound;
+                  numberChanged++;
+                  // check infeasible (relaxed)
+                  if (nowUpper - newBound < 
+                      -100.0*tolerance) {
+                    numberInfeasible++;
+                  }
+                  // adjust
+                  double now;
+                  if (nowLower<-large) {
+                    now=0.0;
+                    infiniteUpper--;
+                  } else {
+                    now = nowLower;
+                  }
+                  maximumUp += (newBound-now) * value;
+                  nowLower = newBound;