Changeset 1676 for trunk

Timestamp:

Jan 20, 2011 12:00:28 PM (9 years ago)

Author:

forrest

Message:

For GUB

Location:

trunk/Clp/src

Files:

• ClpDynamicExampleMatrix.cpp (modified) (1 diff)
• ClpDynamicMatrix.cpp (modified) (30 diffs)
• ClpDynamicMatrix.hpp (modified) (2 diffs)
• ClpNonLinearCost.cpp (modified) (3 diffs)
• ClpSimplex.cpp (modified) (4 diffs)
• ClpSimplexOther.cpp (modified) (2 diffs)
• ClpSimplexOther.hpp (modified) (1 diff)
• ClpSimplexPrimal.cpp (modified) (3 diffs)
• ClpSolve.cpp (modified) (1 diff)

trunk/Clp/src/ClpDynamicExampleMatrix.cpp

@@ -655 +655 @@
           }
      }
-     ClpDynamicMatrix::createVariable(model, bestSequence);
+     ClpDynamicMatrix::createVariable(model, bestSequence/*, bestSequence2*/);
      // clear for next iteration
      savedBestSequence_ = -1;

trunk/Clp/src/ClpDynamicMatrix.cpp

@@ -94 +94 @@
      lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
      upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
-     status_ = ClpCopyOfArray(rhs.status_, numberSets_);
+     status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
      model_ = rhs.model_;
      sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
@@ -110 +110 @@
      maximumGubColumns_ = rhs.maximumGubColumns_;
      maximumElements_ = rhs.maximumElements_;
-     startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_);
+     startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
      next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
      startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
@@ -119 +119 @@
      columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
      columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
-     dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);
+     dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
 }
 
@@ -143 +143 @@
      else
           maximumElements_ = 0;
-     startSet_ = new int [numberSets_];
+     startSet_ = new int [numberSets_+1];
      next_ = new int [maximumGubColumns_];
      // fill in startSet and next
@@ -157 +157 @@
           }
      }
+     startSet_[numberSets_] = starts[numberSets_];
      int numberColumns = model->numberColumns();
      int numberRows = model->numberRows();
@@ -165 +166 @@
      int frequency = model->factorizationFrequency();
      int numberGubInSmall = numberRows + frequency + CoinMin(frequency, numberSets_) + 4;
+     // But we may have two per row + one for incoming (make it two) 
+     numberGubInSmall = CoinMax(2*numberRows+2,numberGubInSmall);
      // for small problems this could be too big
      //numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
@@ -256 +259 @@
      keyVariable_ = new int[numberSets_];
      if (status) {
-          status_ = ClpCopyOfArray(status, numberSets_);
+          status_ = ClpCopyOfArray(status, 2*numberSets_);
           assert (dynamicStatus);
-          dynamicStatus_ = ClpCopyOfArray(dynamicStatus, numberGubColumns_);
+          dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
      } else {
-          status_ = new unsigned char [numberSets_];
+          status_ = new unsigned char [2*numberSets_];
           memset(status_, 0, numberSets_);
           int i;
@@ -267 +270 @@
                setStatus(i, ClpSimplex::basic);
           }
-          dynamicStatus_ = new unsigned char [numberGubColumns_];
+          dynamicStatus_ = new unsigned char [2*numberGubColumns_];
           memset(dynamicStatus_, 0, numberGubColumns_); // for clarity
           for (i = 0; i < numberGubColumns_; i++)
@@ -353 +356 @@
           lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
           upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
-          status_ = ClpCopyOfArray(rhs.status_, numberSets_);
+          status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
           model_ = rhs.model_;
           sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
@@ -369 +372 @@
           maximumGubColumns_ = rhs.maximumGubColumns_;
           maximumElements_ = rhs.maximumElements_;
-          startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_);
+          startSet_ = ClpCopyOfArray(rhs.startSet_, numberSets_+1);
           next_ = ClpCopyOfArray(rhs.next_, maximumGubColumns_);
           startColumn_ = ClpCopyOfArray(rhs.startColumn_, maximumGubColumns_ + 1);
@@ -378 +381 @@
           columnLower_ = ClpCopyOfArray(rhs.columnLower_, maximumGubColumns_);
           columnUpper_ = ClpCopyOfArray(rhs.columnUpper_, maximumGubColumns_);
-          dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, maximumGubColumns_);
+          dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_, 2*maximumGubColumns_);
      }
      return *this;
@@ -533 +536 @@
                                    // check flagged variable and correct dj
                                    if (!flagged(iSequence)) {
+                                     if (false/*status == atLowerBound
+                                                &&keyValue(iSet)<1.0e-7*/) {
+                                       // can't come in because
+                                       // of ones at ub
+                                       numberWanted++;
+                                     } else {
+                                     
                                         bestDj = value;
                                         bestSequence = structuralOffset + iSequence;
                                         bestDjMod = djMod;
                                         bestSet = iSet;
+                                     }
                                    } else {
                                         // just to make sure we don't exit before got something
@@ -576 +587 @@
                            )
 {
-     //forceRefresh=true;
+     // forceRefresh=true;printf("take out forceRefresh\n");
      if (!model_->numberIterations())
           forceRefresh = true;
@@ -735 +746 @@
                                    } else if (getDynamicStatus(j) == atUpperBound) {
                                         value = columnUpper_[j];
+                                        assert (value<1.0e30);
                                    } else if (getDynamicStatus(j) == soloKey) {
                                         value = keyValue(iSet);
-                                   }
+                                   }
                                    objectiveOffset += value * cost_[j];
                               }
@@ -753 +765 @@
                     for (iSet = 0; iSet < numberSets_; iSet++) {
                          int kRow = toIndex_[iSet];
-                         if (kRow >= 0)
-                              kRow += numberStaticRows_;
-                         int j = startSet_[iSet];
-                         while (j >= 0) {
-                              double value = solution[j];
-                              if (value) {
-                                   for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
-                                        int iRow = row_[k];
-                                        rhsOffset_[iRow] -= element_[k] * value;
-                                   }
-                                   if (kRow >= 0)
-                                        rhsOffset_[kRow] -= value;
-                              }
-                              j = next_[j]; //onto next in set
-                         }
+                         if (kRow >= 0) 
+                           kRow += numberStaticRows_;
+                         int j = startSet_[iSet];
+                         while (j >= 0) {
+                           //? use DynamicStatus status = getDynamicStatus(j);
+                           double value = solution[j];
+                           if (value) {
+                             for (CoinBigIndex k = startColumn_[j]; k < startColumn_[j+1]; k++) {
+                               int iRow = row_[k];
+                               rhsOffset_[iRow] -= element_[k] * value;
+                             }
+                             if (kRow >= 0)
+                               rhsOffset_[kRow] -= value;
+                           }
+                           j = next_[j]; //onto next in set
+                         }
                     }
                     delete [] solution;
@@ -1073 +1086 @@
      // save status
      case 5: {
+       memcpy(status_+numberSets_,status_,numberSets_);
+       memcpy(dynamicStatus_+maximumGubColumns_,
+              dynamicStatus_,maximumGubColumns_);
      }
      break;
      // restore status
      case 6: {
-          // maybe mismatch - redo all in small model
-          //for (int i=firstDynamic_;i<firstAvailable_;i++) {
-          //int jColumn = id_[i-firstDynamic_];
-          //setDynamicStatus(jColumn,inSmall);
-          //}
+       memcpy(status_,status_+numberSets_,numberSets_);
+       memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
+              maximumGubColumns_);
+       initialProblem();
      }
      break;
@@ -1172 +1187 @@
           // first flag
           if (number >= firstDynamic_ && number < lastDynamic_) {
-               assert (number == model->sequenceIn());
-               // id will be sitting at firstAvailable
-               int sequence = id_[firstAvailable_-firstDynamic_];
-               setFlagged(sequence);
-               model->clearFlagged(firstAvailable_);
-          }
+            int sequence = id_[number-firstDynamic_];
+            setFlagged(sequence);
+            //model->clearFlagged(number);
+          } else if (number>=model_->numberColumns()+numberStaticRows_) {
+            // slack
+            int iSet = fromIndex_[number-model_->numberColumns()-
+                                  numberStaticRows_];
+            setFlaggedSlack(iSet);
+            //model->clearFlagged(number);
+          }
      }
      case 11: {
-          int sequenceIn = model->sequenceIn();
-          if (sequenceIn >= firstDynamic_ && sequenceIn < lastDynamic_) {
-               assert (number == model->sequenceIn());
+          //int sequenceIn = model->sequenceIn();
+          if (number >= firstDynamic_ && number < lastDynamic_) {
+            //assert (number == model->sequenceIn());
                // take out variable (but leave key)
                double * cost = model->costRegion();
@@ -1190 +1209 @@
                int * length = matrix_->getMutableVectorLengths();
 #ifndef NDEBUG
-               if (length[sequenceIn]) {
+               if (length[number]) {
                     int * row = matrix_->getMutableIndices();
                     CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
-                    int iRow = row[startColumn[sequenceIn] + length[sequenceIn] - 1];
+                    int iRow = row[startColumn[number] + length[number] - 1];
                     int which = iRow - numberStaticRows_;
                     assert (which >= 0);
@@ -1210 +1229 @@
 
                // not really in small problem
-               int iBig = id_[sequenceIn-firstDynamic_];
-               if (model->getStatus(sequenceIn) == ClpSimplex::atLowerBound) {
+               int iBig = id_[number-firstDynamic_];
+               if (model->getStatus(number) == ClpSimplex::atLowerBound) {
                     setDynamicStatus(iBig, atLowerBound);
                     if (columnLower_)
-                         modifyOffset(sequenceIn, columnLower_[iBig]);
+                         modifyOffset(number, columnLower_[iBig]);
                } else {
                     setDynamicStatus(iBig, atUpperBound);
-                    modifyOffset(sequenceIn, columnUpper_[iBig]);
-               }
+                    modifyOffset(number, columnUpper_[iBig]);
+               }
+          } else if (number>=model_->numberColumns()+numberStaticRows_) {
+            // slack
+            int iSet = fromIndex_[number-model_->numberColumns()-
+                                  numberStaticRows_];
+            printf("what now - set %d\n",iSet);
           }
      }
@@ -1302 +1326 @@
                double value = solution[iColumn];
                bool toLowerBound = true;
+               assert (jColumn>=0);assert (iColumn>=0);
                if (columnUpper_) {
                     if (!columnLower_) {
                          if (fabs(value - columnUpper_[jColumn]) < fabs(value))
                               toLowerBound = false;
-                    } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[iColumn])) {
+                    } else if (fabs(value - columnUpper_[jColumn]) < fabs(value - columnLower_[jColumn])) {
                          toLowerBound = false;
                     }
@@ -1509 +1534 @@
                n++;
           }
+          int k=0;
+          for (int j=startSet_[i];j<startSet_[i+1];j++) {
+            if (getDynamicStatus(j)==soloKey)
+              k++;
+          }
+          assert (k<2);
      }
      assert (n == numberSets_);
@@ -1597 +1628 @@
                     j = next_[j]; //onto next in set
                }
-          }
+#if 0
+               // slack must be feasible
+               double oldValue=value;
+               value = CoinMax(value,lowerSet_[iSet]);
+               value = CoinMin(value,upperSet_[iSet]);
+               if (value!=oldValue)
+                 printf("using %g (not %g) for slack on set %d (%g,%g)\n",
+                        value,oldValue,iSet,lowerSet_[iSet],upperSet_[iSet]);
+#endif
+            }
      }
      return value;
@@ -1676 +1716 @@
                     model->setStatus(iSequence, getStatus(savedBestSet_));
                     model->djRegion()[iSequence] = savedBestGubDual_;
-                    if (getStatus(savedBestSet_) == ClpSimplex::atUpperBound)
-                         solution[iSequence] = columnUpper[iSequence];
-                    else
-                         solution[iSequence] = columnLower[iSequence];
-                    // create variable and pivot in
+                    solution[iSequence] = valueOfKey;
+                    // create variable and pivot in
                     int key = keyVariable_[savedBestSet_];
                     setDynamicStatus(key, inSmall);
@@ -1770 +1807 @@
                     model->setStatus(iSequence, ClpSimplex::basic);
                     model->djRegion()[iSequence] = 0.0;
-                    solution[iSequence] = shift;
+                    solution[iSequence] = valueOfKey+shift;
                     rhsOffset_[newRow] = -shift; // sign?
                }
@@ -1816 +1853 @@
                //printf("best %d\n",bestSequence2);
                model->solutionRegion()[firstAvailable_] = 0.0;
+               model->clearFlagged(firstAvailable_);
                if (!columnLower_ && !columnUpper_) {
                     model->setStatus(firstAvailable_, ClpSimplex::atLowerBound);
@@ -2109 +2147 @@
           int numberActive = 0;
           int whichKey = -1;
-          if (getStatus(iSet) == ClpSimplex::basic)
+          if (getStatus(iSet) == ClpSimplex::basic) {
                whichKey = maximumGubColumns_ + iSet;
-          else
+               numberActive = 1;
+          } else {
                whichKey = -1;
+          }
           int j = startSet_[iSet];
           while (j >= 0) {
                assert (getDynamicStatus(j) != soloKey || whichKey == -1);
-               if (getDynamicStatus(j) == inSmall)
+               if (getDynamicStatus(j) == inSmall) {
                     numberActive++;
-               else if (getDynamicStatus(j) == soloKey)
+               } else if (getDynamicStatus(j) == soloKey) {
                     whichKey = j;
+                    numberActive++;
+               }
                j = next_[j]; //onto next in set
           }
-          if (getStatus(iSet) == ClpSimplex::basic && numberActive)
-               numberActive++;
-          if (numberActive) {
-               assert (numberActive > 1);
+          if (numberActive > 1) {
                int iRow = numberActiveSets_ + numberStaticRows_;
                rowSolution[iRow] = 0.0;
@@ -2189 +2228 @@
                               else
                                    columnUpper[firstAvailable_] = COIN_DBL_MAX;
-                              if (status == atLowerBound) {
+                              if (status != atUpperBound) {
                                    solution[firstAvailable_] = columnLower[firstAvailable_];
                               } else {
@@ -2203 +2242 @@
                toIndex_[iSet] = numberActiveSets_;
                fromIndex_[numberActiveSets_++] = iSet;
+          } else {
+            // solo key
+            bool needKey;
+            if (numberActive) {
+              if (whichKey<maximumGubColumns_) {
+                // structural - assume ok
+                needKey = false;
+              } else {
+                // slack
+                keyVariable_[iSet] = maximumGubColumns_ + iSet;
+                double value = keyValue(iSet);
+                if (value<lowerSet_[iSet]-1.0e-8||
+                    value>upperSet_[iSet]+1.0e-8)
+                  needKey=true;
+              }
+            } else {
+              needKey = true;
+            }
+            if (needKey) {
+              // all to lb then up some (slack/null if possible)
+              int length=99999999;
+              int which=-1;
+              double sum=0.0;
+              for (int iColumn=startSet_[iSet];iColumn<startSet_[iSet+1];iColumn++) {
+                setDynamicStatus(iColumn,atLowerBound);
+                sum += columnLower_[iColumn];
+                if (length>startColumn_[iColumn+1]-startColumn_[iColumn]) {
+                  which=iColumn;
+                  length=startColumn_[iColumn+1]-startColumn_[iColumn];
+                }
+              }
+              if (sum>lowerSet_[iSet]-1.0e-8) {
+                // slack can be basic
+                setStatus(iSet,ClpSimplex::basic);
+                keyVariable_[iSet] = maximumGubColumns_ + iSet;
+              } else {
+                // use shortest
+                setDynamicStatus(which,soloKey);
+                keyVariable_[iSet] = which;
+                setStatus(iSet,ClpSimplex::atLowerBound);
+              }
+            }
           }
           assert (toIndex_[iSet] >= 0 || whichKey >= 0);
           keyVariable_[iSet] = whichKey;
      }
+     numberActiveColumns_ = firstAvailable_;
      return;
+}
+// Writes out model (without names)
+void 
+ClpDynamicMatrix::writeMps(const char * name)
+{
+  int numberTotalRows = numberStaticRows_+numberSets_;
+  int numberTotalColumns = firstDynamic_+numberGubColumns_;
+  // over estimate
+  int numberElements = getNumElements()+startColumn_[numberGubColumns_]
+    + numberGubColumns_;
+  double * columnLower = new double [numberTotalColumns];
+  double * columnUpper = new double [numberTotalColumns];
+  double * cost = new double [numberTotalColumns];
+  double * rowLower = new double [numberTotalRows];
+  double * rowUpper = new double [numberTotalRows];
+  CoinBigIndex * start = new CoinBigIndex[numberTotalColumns+1];
+  int * row = new int [numberElements];
+  double * element = new double [numberElements];
+  // Fill in
+  const CoinBigIndex * startA = getVectorStarts();
+  const int * lengthA = getVectorLengths();
+  const int * rowA = getIndices();
+  const double * elementA = getElements();
+  const double * columnLowerA = model_->columnLower();
+  const double * columnUpperA = model_->columnUpper();
+  const double * costA = model_->objective();
+  const double * rowLowerA = model_->rowLower();
+  const double * rowUpperA = model_->rowUpper();
+  start[0]=0;
+  numberElements=0;
+  for (int i=0;i<firstDynamic_;i++) {
+    columnLower[i] = columnLowerA[i];
+    columnUpper[i] = columnUpperA[i];
+    cost[i] = costA[i];
+    for (CoinBigIndex j = startA[i];j<startA[i]+lengthA[i];j++) {
+      row[numberElements] = rowA[j];
+      element[numberElements++]=elementA[j];
+    }
+    start[i+1]=numberElements;
+  }
+  for (int i=0;i<numberStaticRows_;i++) {
+    rowLower[i] = rowLowerA[i];
+    rowUpper[i] = rowUpperA[i];
+  }
+  int putC=firstDynamic_;
+  int putR=numberStaticRows_;
+  for (int i=0;i<numberSets_;i++) {
+    rowLower[putR]=lowerSet_[i];
+    rowUpper[putR]=upperSet_[i];
+    for (CoinBigIndex k=startSet_[i];k<startSet_[i+1];k++) {
+      columnLower[putC]=columnLower_[k];
+      columnUpper[putC]=columnUpper_[k];
+      cost[putC]=cost_[k];
+      putC++;
+      for (CoinBigIndex j = startColumn_[k];j<startColumn_[k+1];j++) {
+        row[numberElements] = row_[j];
+        element[numberElements++]=element_[j];
+      }
+      row[numberElements] = putR;
+      element[numberElements++]=1.0;
+      start[putC]=numberElements;
+    }
+    putR++;
+  }
+
+  assert (putR==numberTotalRows);
+  assert (putC==numberTotalColumns);
+  ClpSimplex modelOut;
+  modelOut.loadProblem(numberTotalColumns,numberTotalRows,
+                       start,row,element,
+                       columnLower,columnUpper,cost,
+                       rowLower,rowUpper);
+  modelOut.writeMps(name);
+  delete [] columnLower;
+  delete [] columnUpper;
+  delete [] cost;
+  delete [] rowLower;
+  delete [] rowUpper;
+  delete [] start;
+  delete [] row; 
+  delete [] element;
 }
 // Adds in a column to gub structure (called from descendant)
@@ -2236 +2399 @@
                     if (columnUpper_ && upper != columnUpper_[j])
                          odd = true;
-                    if (odd)
+                    if (odd) {
                          printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
                                 cost, lower, upper, cost_[j],
                                 columnLower_ ? columnLower_[j] : 0.0,
                                 columnUpper_ ? columnUpper_[j] : 1.0e100);
-                    else
+                    } else {
+                         setDynamicStatus(j, status);
                          return j;
+                    }
                }
           }
@@ -2259 +2424 @@
           for (i = 0; i < numberGubColumns_; i++) {
                CoinBigIndex end = startColumn_[i+1];
+               // what about ubs if column generation?
                if (getDynamicStatus(i) != atLowerBound) {
                     // keep in

trunk/Clp/src/ClpDynamicMatrix.hpp

@@ -95 +95 @@
      /// Does gub crash
      void gubCrash();
+     /// Writes out model (without names)
+     void writeMps(const char * name);
      /// Populates initial matrix from dynamic status
      void initialProblem();
@@ -166 +168 @@
           st_byte = static_cast<unsigned char>(st_byte | status);
      }
+     /// Whether flagged slack
+     inline bool flaggedSlack(int i) const {
+          return (status_[i] & 8) != 0;
+     }
+     inline void setFlaggedSlack(int i) {
+          status_[i] = static_cast<unsigned char>(status_[i] | 8);
+     }
+     inline void unsetFlaggedSlack(int i) {
+          status_[i] = static_cast<unsigned char>(status_[i] & ~8);
+     }
      /// Number of sets (dynamic rows)
      inline int numberSets() const {
           return numberSets_;
      }
+     /// Number of possible gub variables
+     inline int numberGubEntries() const
+     { return startSet_[numberSets_];}
+     /// Sets
+     inline int * startSets() const
+     { return startSet_;}
      /// Whether flagged
      inline bool flagged(int i) const {

trunk/Clp/src/ClpNonLinearCost.cpp

@@ -608 +608 @@
                     }
                }
+               //#define PRINT_DETAIL7 2
+#if PRINT_DETAIL7>1
+               printf("NL %d sol %g bounds %g %g\n",
+                      iSequence,solution[iSequence],
+                      lowerValue,upperValue);
+#endif
                switch(status) {
 
@@ -625 +631 @@
                                         printf("nonlincostb %d %g %g %g\n",
                                                iSequence, lowerValue, solution[iSequence], lower_[iRange+2]);
+#endif
+#if PRINT_DETAIL7
+                                   printf("**NL %d sol %g below %g\n",
+                                          iSequence,solution[iSequence],
+                                          lowerValue);
 #endif
                                    sumInfeasibilities_ += value;
@@ -642 +653 @@
                                         printf("nonlincostu %d %g %g %g\n",
                                                iSequence, lower_[iRange-1], solution[iSequence], upperValue);
+#endif
+#if PRINT_DETAIL7
+                                   printf("**NL %d sol %g above %g\n",
+                                          iSequence,solution[iSequence],
+                                          upperValue);
 #endif
                                    sumInfeasibilities_ += value;

trunk/Clp/src/ClpSimplex.cpp

@@ -1933 +1933 @@
      int cycle = progress_.cycle(in, out,
                                  directionIn_, directionOut_);
-     if (cycle > 0 && objective_->type() < 2) {
+     if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
           //if (cycle>0) {
           if (handler_->logLevel() >= 63)
@@ -1945 +1945 @@
                forceFactorization_ = CoinMax(1, cycle - off[extra]);
           } else {
-               // need to reject something
+            /* need to reject something
+               should be better if don't reject incoming
+               as it is in basis */
                int iSequence;
-               if (algorithm_ > 0)
-                    iSequence = sequenceIn_;
-               else
+               //if (algorithm_ > 0)
+               //   iSequence = sequenceIn_;
+               //else
                     iSequence = sequenceOut_;
                char x = isColumn(iSequence) ? 'C' : 'R';
@@ -1993 +1995 @@
                forceFactorization_ = -1; //off
           return 1;
-     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
+     } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
           double random = randomNumberGenerator_.randomDouble();
           int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
@@ -8250 +8252 @@
           whatsChanged_ &= ~0xffff;
      }
+     double saveObjValue = objectiveValue_;
      deleteRim(getRidOfData);
+     if (matrix_->type()>=15)
+       objectiveValue_ = saveObjValue;
      // Skip message if changing algorithms
      if (problemStatus_ != 10) {

trunk/Clp/src/ClpSimplexOther.cpp

@@ -16 +16 @@
 #include "ClpFactorization.hpp"
 #include "ClpDualRowDantzig.hpp"
+#include "ClpDynamicMatrix.hpp"
 #include "CoinPackedMatrix.hpp"
 #include "CoinIndexedVector.hpp"
@@ -4068 +4069 @@
      checkSolutionInternal();
 }
+// Returns gub version of model or NULL
+ClpSimplex * 
+ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns,
+                            int neededGub,
+                            int factorizationFrequency)
+{
+  // find gub
+  int numberRows = this->numberRows();
+  int numberColumns = this->numberColumns();
+  int iRow, iColumn;
+  int * columnIsGub = new int [numberColumns];
+  const double * columnLower = this->columnLower();
+  const double * columnUpper = this->columnUpper();
+  int numberFixed=0;
+  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
+    if (columnUpper[iColumn] == columnLower[iColumn]) {
+      columnIsGub[iColumn]=-2;
+      numberFixed++;
+    } else if (columnLower[iColumn]>=0) {
+      columnIsGub[iColumn]=-1;
+    } else {
+      columnIsGub[iColumn]=-3;
+    }
+  }
+  CoinPackedMatrix * matrix = this->matrix();
+  // get row copy
+  CoinPackedMatrix rowCopy = *matrix;
+  rowCopy.reverseOrdering();
+  const int * column = rowCopy.getIndices();
+  const int * rowLength = rowCopy.getVectorLengths();
+  const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
+  const double * element = rowCopy.getElements();
+  int numberNonGub = 0;
+  int numberEmpty = numberRows;
+  int * rowIsGub = new int [numberRows];
+  int smallestGubRow=-1;
+  int count=numberColumns+1;
+  double * rowLower = this->rowLower();
+  double * rowUpper = this->rowUpper();
+  // make sure we can get rid of upper bounds
+  double * fixedRow = new double [numberRows];
+  for (iRow = 0 ; iRow < numberRows ; iRow++) {
+    double sumFixed=0.0;
+    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
+      int iColumn = column[j];
+      double value = columnLower[iColumn];
+      if (value) 
+        sumFixed += element[j] * value;
+    }
+    fixedRow[iRow]=rowUpper[iRow]-sumFixed;
+  }
+  for (iRow = numberRows - 1; iRow >= 0; iRow--) {
+    bool gubRow = true;
+    int numberInRow=0;
+    double sumFixed=0.0;
+    double gap = fixedRow[iRow]-1.0e-12;
+    for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
+      int iColumn = column[j];
+      if (columnIsGub[iColumn]!=-2) {
+        if (element[j] != 1.0||columnIsGub[iColumn]==-3||
+            columnUpper[iColumn]-columnLower[iColumn]<gap) {
+          gubRow = false;
+          break;
+        } else {
+          numberInRow++;
+          if (columnIsGub[iColumn] >= 0) {
+            gubRow = false;
+            break;
+          }
+        }
+      } else {
+        sumFixed += columnLower[iColumn]*element[j];
+      }
+    }
+    if (!gubRow) {
+      whichRows[numberNonGub++] = iRow;
+      rowIsGub[iRow] = -1;
+    } else if (numberInRow) {
+      if (numberInRow<count) {
+        count = numberInRow;
+        smallestGubRow=iRow;
+      }
+      for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
+        int iColumn = column[j];
+        if (columnIsGub[iColumn]!=-2) 
+          columnIsGub[iColumn] = iRow;
+      }
+      rowIsGub[iRow] = 0;
+    } else {
+      // empty row!
+      whichRows[--numberEmpty] = iRow;
+      rowIsGub[iRow] = -2;
+      if (sumFixed>rowUpper[iRow]+1.0e-4||
+          sumFixed<rowLower[iRow]-1.0e-4) {
+        fprintf(stderr,"******** No infeasible empty rows - please!\n");
+        abort();
+      }
+    }
+  }
+  delete [] fixedRow;
+  char message[100];
+  int numberGub = numberEmpty - numberNonGub;
+  if (numberGub >= neededGub) { 
+    sprintf(message,"%d gub rows", numberGub);
+    handler_->message(CLP_GENERAL, messages_)
+      << message << CoinMessageEol;
+    int numberNormal = 0;
+    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
+      if (columnIsGub[iColumn] < 0  && columnIsGub[iColumn] !=-2) {
+        whichColumns[numberNormal++] = iColumn;
+      }
+    }
+    if (!numberNormal) {
+      sprintf(message,"Putting back one gub row to make non-empty");
+      handler_->message(CLP_GENERAL, messages_)
+        << message << CoinMessageEol;
+      rowIsGub[smallestGubRow]=-1;
+      whichRows[numberNonGub++] = smallestGubRow;
+      for (int j = rowStart[smallestGubRow]; 
+           j < rowStart[smallestGubRow] + rowLength[smallestGubRow]; j++) {
+        int iColumn = column[j];
+        if (columnIsGub[iColumn]>=0) {
+          columnIsGub[iColumn]=-4;
+          whichColumns[numberNormal++] = iColumn;
+        }
+      }
+    }
+    std::sort(whichRows,whichRows+numberNonGub);
+    std::sort(whichColumns,whichColumns+numberNormal);
+    double * lower = CoinCopyOfArray(this->rowLower(),numberRows);
+    double * upper = CoinCopyOfArray(this->rowUpper(),numberRows);
+    // leave empty rows at end
+    numberEmpty = numberRows-numberEmpty;
+    const int * row = matrix->getIndices();
+    const int * columnLength = matrix->getVectorLengths();
+    const CoinBigIndex * columnStart = matrix->getVectorStarts();
+    const double * elementByColumn = matrix->getElements();
+    // Fixed at end
+    int put2 = numberColumns-numberFixed;
+    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
+      if (columnIsGub[iColumn] ==-2) {
+        whichColumns[put2++] = iColumn;
+        double value = columnLower[iColumn];
+        for (int j = columnStart[iColumn]; 
+             j < columnStart[iColumn] + columnLength[iColumn]; j++) {
+          int iRow = row[j];
+          if (lower[iRow]>-1.0e20)
+            lower[iRow] -= value*element[j];
+          if (upper[iRow]<1.0e20)
+            upper[iRow] -= value*element[j];
+        }
+      }
+    }
+    int put = numberNormal;
+    ClpSimplex * model2 =
+      new ClpSimplex(this, numberNonGub, whichRows , numberNormal, whichColumns);
+    // scale
+    double * scaleArray = new double [numberRows];
+    for (int i=0;i<numberRows;i++) {
+      scaleArray[i]=1.0;
+      if (rowIsGub[i]==-1) {
+        double largest = 1.0e-30;
+        double smallest = 1.0e30;
+        for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
+          int iColumn = column[j];
+          if (columnIsGub[iColumn]!=-2) {
+            double value =fabs(element[j]);
+            largest = CoinMax(value,largest);
+            smallest = CoinMin(value,smallest);
+          }
+        }
+        double scale = CoinMax(0.001,1.0/sqrt(largest*smallest));
+        scaleArray[i]=scale;
+        if (lower[i]>-1.0e30)
+          lower[i] *= scale;
+        if (upper[i]<1.0e30)
+          upper[i] *= scale;
+      }
+    }
+    // scale partial matrix
+    {
+      CoinPackedMatrix * matrix = model2->matrix();
+      const int * row = matrix->getIndices();
+      const int * columnLength = matrix->getVectorLengths();
+      const CoinBigIndex * columnStart = matrix->getVectorStarts();
+      double * element = matrix->getMutableElements();
+      for (int i=0;i<numberNormal;i++) {
+        for (int j = columnStart[i]; 
+             j < columnStart[i] + columnLength[i]; j++) {
+          int iRow = row[j];
+          iRow = whichRows[iRow];
+          double scaleBy = scaleArray[iRow];
+          element[j] *= scaleBy;
+        }
+      }
+    }
+    // adjust rhs
+    double * rowLower = model2->rowLower();
+    double * rowUpper = model2->rowUpper();
+    for (int i=0;i<numberNonGub;i++) {
+      int iRow = whichRows[i];
+      rowLower[i] = lower[iRow];
+      rowUpper[i] = upper[iRow];
+    }
+    int numberGubColumns = numberColumns - put - numberFixed;
+    CoinBigIndex numberElements=0;
+    int * temp1 = new int [numberRows+1];
+    // get counts
+    memset(temp1,0,numberRows*sizeof(int));
+    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
+      int iGub = columnIsGub[iColumn];
+      if (iGub>=0) {
+        numberElements += columnLength[iColumn]-1;
+        temp1[iGub]++;
+      }
+    }
+    /* Optional but means can eventually simplify coding
+       could even add in fixed slacks to deal with
+       singularities - but should not be necessary */
+    int numberSlacks=0;
+    for (int i = 0; i < numberRows; i++) {
+      if (rowIsGub[i]>=0) {
+        if (lower[i]<upper[i]) {
+          numberSlacks++;
+          temp1[i]++;
+        }
+      }
+    }
+    int * gubStart = new int [numberGub+1];
+    numberGub=0;
+    gubStart[0]=0;
+    for (int i = 0; i < numberRows; i++) {
+      if (rowIsGub[i]>=0) {
+        rowIsGub[i]=numberGub;
+        gubStart[numberGub+1]=gubStart[numberGub]+temp1[i];
+        temp1[numberGub]=0;
+        lower[numberGub]=lower[i];
+        upper[numberGub]=upper[i];
+        whichRows[numberNonGub+numberGub]=i;
+        numberGub++;
+      }
+    }
+    int numberGubColumnsPlus = numberGubColumns + numberSlacks;
+    double * lowerColumn2 = new double [numberGubColumnsPlus];
+    CoinFillN(lowerColumn2, numberGubColumnsPlus, 0.0);
+    double * upperColumn2 = new double [numberGubColumnsPlus];
+    CoinFillN(upperColumn2, numberGubColumnsPlus, COIN_DBL_MAX);
+    int * start2 = new int[numberGubColumnsPlus+1];
+    int * row2 = new int[numberElements];
+    double * element2 = new double[numberElements];
+    double * cost2 = new double [numberGubColumnsPlus];
+    CoinFillN(cost2, numberGubColumnsPlus, 0.0);
+    const double * cost = this->objective();
+    put = numberNormal;
+    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
+      int iGub = columnIsGub[iColumn];
+      if (iGub>=0) {
+        // TEMP
+        //this->setColUpper(iColumn,COIN_DBL_MAX);
+        iGub = rowIsGub[iGub];
+        assert (iGub>=0);
+        int kPut = put+gubStart[iGub]+temp1[iGub];
+        temp1[iGub]++;
+        whichColumns[kPut]=iColumn;
+      }
+    }
+    for (int i = 0; i < numberRows; i++) {
+      if (rowIsGub[i]>=0) {
+        int iGub = rowIsGub[i];
+        if (lower[iGub]<upper[iGub]) {
+          int kPut = put+gubStart[iGub]+temp1[iGub];
+          temp1[iGub]++;
+          whichColumns[kPut]=iGub+numberColumns;
+        }
+      }
+    }
+    //this->primal(1); // TEMP
+    // redo rowIsGub to give lookup
+    for (int i=0;i<numberRows;i++)
+      rowIsGub[i]=-1;
+    for (int i=0;i<numberNonGub;i++) 
+      rowIsGub[whichRows[i]]=i;
+    start2[0]=0;
+    numberElements = 0;
+    for (int i=0;i<numberGubColumnsPlus;i++) {
+      int iColumn = whichColumns[put++];
+      if (iColumn<numberColumns) {
+        cost2[i] = cost[iColumn];
+        lowerColumn2[i] = columnLower[iColumn];
+        upperColumn2[i] = columnUpper[iColumn];
+        upperColumn2[i] = COIN_DBL_MAX; 
+        for (int j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
+          int iRow = row[j];
+          double scaleBy = scaleArray[iRow];
+          iRow = rowIsGub[iRow];
+          if (iRow >= 0) {
+            row2[numberElements] = iRow;
+            element2[numberElements++] = elementByColumn[j]*scaleBy;
+          }
+        }
+      } else {
+        // slack
+        int iGub = iColumn-numberColumns;
+        double slack = upper[iGub]-lower[iGub];
+        assert (upper[iGub]<1.0e20);
+        lower[iGub]=upper[iGub];
+        cost2[i] = 0;
+        lowerColumn2[i] = 0;
+        upperColumn2[i] = slack;
+        upperColumn2[i] = COIN_DBL_MAX;
+      }
+      start2[i+1] = numberElements;
+    }
+    // clean up bounds on variables
+    for (int iSet=0;iSet<numberGub;iSet++) {
+      double lowerValue=0.0;
+      for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
+        lowerValue += lowerColumn2[i];
+      }
+      assert (lowerValue<upper[iSet]+1.0e-6);
+      double gap = CoinMax(0.0,upper[iSet]-lowerValue);
+      for (int i=gubStart[iSet];i<gubStart[iSet+1];i++) {
+        if (upperColumn2[i]<1.0e30) {
+          upperColumn2[i] = CoinMin(upperColumn2[i],
+                                    lowerColumn2[i]+gap);
+        }
+      }
+    }
+    sprintf(message,"** Before adding matrix there are %d rows and %d columns",
+           model2->numberRows(), model2->numberColumns());
+    handler_->message(CLP_GENERAL, messages_)
+      << message << CoinMessageEol;
+    delete [] scaleArray;
+    delete [] temp1;
+    model2->setFactorizationFrequency(factorizationFrequency);
+    ClpDynamicMatrix * newMatrix = 
+      new ClpDynamicMatrix(model2, numberGub,
+                                  numberGubColumnsPlus, gubStart,
+                                  lower, upper,
+                                  start2, row2, element2, cost2,
+                                  lowerColumn2, upperColumn2);
+    delete [] gubStart;
+    delete [] lowerColumn2;
+    delete [] upperColumn2;
+    delete [] start2;
+    delete [] row2;
+    delete [] element2;
+    delete [] cost2;
+    delete [] lower;
+    delete [] upper;
+    model2->replaceMatrix(newMatrix,true);
+#ifdef EVERY_ITERATION
+    {
+      ClpDynamicMatrix * gubMatrix =
+        dynamic_cast< ClpDynamicMatrix*>(model2->clpMatrix());
+      assert(gubMatrix);
+      gubMatrix->writeMps("gub.mps");
+    }
+#endif
+    delete [] columnIsGub;
+    delete [] rowIsGub;
+    newMatrix->switchOffCheck();
+#ifdef EVERY_ITERATION
+    newMatrix->setRefreshFrequency(1/*000*/);
+#else
+    newMatrix->setRefreshFrequency(1000);
+#endif
+    sprintf(message,
+            "** While after adding matrix there are %d rows and %d columns",
+            model2->numberRows(), model2->numberColumns());
+    handler_->message(CLP_GENERAL, messages_)
+      << message << CoinMessageEol;
+    model2->setSpecialOptions(4);    // exactly to bound
+    // Scaling off (done by hand)
+    model2->scaling(0);
+    return model2;
+  } else {
+    delete [] columnIsGub;
+    delete [] rowIsGub;
+    return NULL;
+  }
+}
+// Sets basis from original
+void 
+ClpSimplexOther::setGubBasis(const ClpSimplex &original,const int * whichRows,
+                             const int * whichColumns)
+{
+  ClpDynamicMatrix * gubMatrix =
+    dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
+  assert(gubMatrix);
+  int numberGubColumns = gubMatrix->numberGubColumns();
+  int numberNormal = gubMatrix->firstDynamic();
+  //int lastOdd = gubMatrix->firstAvailable();
+  //int numberTotalColumns = numberNormal + numberGubColumns;
+  //assert (numberTotalColumns==numberColumns+numberSlacks);
+  int numberRows = original.numberRows();
+  int numberColumns = original.numberColumns();
+  int * columnIsGub = new int [numberColumns];
+  int numberNonGub = gubMatrix->numberStaticRows();
+  //assert (firstOdd==numberNormal);
+  double * solution = primalColumnSolution();
+  const double * originalSolution = original.primalColumnSolution();
+  int numberSets = gubMatrix->numberSets();
+  const int * startSet = gubMatrix->startSets();
+  const CoinBigIndex * startColumn = gubMatrix->startColumn();
+  const double * columnLower = gubMatrix->columnLower();
+  for (int i=0;i<numberSets;i++) {
+    for (int j=startSet[i];j<startSet[i+1];j++) {
+      int iColumn = whichColumns[j+numberNormal];
+      if (iColumn<numberColumns)
+        columnIsGub[iColumn] = whichRows[numberNonGub+i];
+    }
+  }
+  int * numberKey = new int [numberRows];
+  memset(numberKey,0,numberRows*sizeof(int));
+  for (int i=0;i<numberGubColumns;i++) {
+    int iOrig = whichColumns[i+numberNormal];
+    if (iOrig<numberColumns) {
+      if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
+        int iRow = columnIsGub[iOrig];
+        assert (iRow>=0);
+        numberKey[iRow]++;
+      }
+    } else {
+      // Set slack
+      int iSet = iOrig - numberColumns;
+      int iRow = whichRows[iSet+numberNonGub];
+      if (original.getRowStatus(iRow)==ClpSimplex::basic) 
+        numberKey[iRow]++;
+    }
+  }
+  /* Before going into cleanMatrix we need
+     gub status set (inSmall just means basic and active)
+     row status set
+  */
+  for (int i = 0; i < numberSets; i++) {
+    gubMatrix->setStatus(i,ClpSimplex::isFixed);
+  }
+  for (int i = 0; i < numberGubColumns; i++) {
+    int iOrig = whichColumns[i+numberNormal];
+    if (iOrig<numberColumns) {
+      ClpSimplex::Status status = original.getColumnStatus(iOrig);
+      if (status==ClpSimplex::atUpperBound) {
+        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atUpperBound);
+      } else if (status==ClpSimplex::atLowerBound) {
+        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
+      } else if (status==ClpSimplex::basic) {
+        int iRow = columnIsGub[iOrig];
+        assert (iRow>=0);
+        assert(numberKey[iRow]);
+        if (numberKey[iRow]==1)
+          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
+        else
+          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
+      }
+    } else {
+      // slack
+      int iSet = iOrig - numberColumns;
+      int iRow = whichRows[iSet+numberNonGub];
+      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
+        assert(numberKey[iRow]);
+        if (numberKey[iRow]==1)
+          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::soloKey);
+        else
+          gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::inSmall);
+      } else {
+        gubMatrix->setDynamicStatus(i,ClpDynamicMatrix::atLowerBound);
+      }
+    }
+  }
+  // deal with sets without key
+  for (int i = 0; i < numberSets; i++) {
+    int iRow = whichRows[numberNonGub+i];
+    if (!numberKey[iRow]) {
+      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
+        gubMatrix->setStatus(i,ClpSimplex::basic);
+      } else { 
+        // If not at lb make key otherwise one with smallest number els
+        double largest=0.0;
+        int fewest=numberRows+1;
+        int chosen=-1;
+        for (int j=startSet[i];j<startSet[i+1];j++) {
+          int length=startColumn[j+1]-startColumn[j];
+          int iOrig = whichColumns[j+numberNormal];
+          double value;
+          if (iOrig<numberColumns) {
+            value=originalSolution[iOrig]-columnLower[j];
+          } else {
+            // slack - take value as 0.0 as will win on length
+            value=0.0;
+          }
+          if (value>largest+1.0e-8) {
+            largest=value;
+            fewest=length;
+            chosen=j;
+          } else if (fabs(value-largest)<=1.0e-8&&length<fewest) {
+            largest=value;
+            fewest=length;
+            chosen=j;
+          }
+        }
+        assert(chosen>=0);
+        // set as key
+        for (int j=startSet[i];j<startSet[i+1];j++) {
+          if (j!=chosen)
+            gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
+          else
+            gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::soloKey);
+        }
+      }
+    }
+  }
+  for (int i = 0; i < numberNormal; i++) {
+    int iOrig = whichColumns[i];
+    setColumnStatus(i,original.getColumnStatus(iOrig));
+    solution[i]=originalSolution[iOrig];
+  }
+  for (int i = 0; i < numberNonGub; i++) {
+    int iOrig = whichRows[i];
+    setRowStatus(i,original.getRowStatus(iOrig));
+  }
+  // Fill in current matrix
+  gubMatrix->initialProblem();
+  delete [] numberKey;
+  delete [] columnIsGub;
+}
+// Restores basis to original
+void 
+ClpSimplexOther::getGubBasis(ClpSimplex &original,const int * whichRows,
+                             const int * whichColumns) const
+{
+  ClpDynamicMatrix * gubMatrix =
+    dynamic_cast< ClpDynamicMatrix*>(clpMatrix());
+  assert(gubMatrix);
+  int numberGubColumns = gubMatrix->numberGubColumns();
+  int numberNormal = gubMatrix->firstDynamic();
+  //int lastOdd = gubMatrix->firstAvailable();
+  //int numberRows = original.numberRows();
+  int numberColumns = original.numberColumns();
+  int numberNonGub = gubMatrix->numberStaticRows();
+  //assert (firstOdd==numberNormal);
+  double * solution = primalColumnSolution();
+  double * originalSolution = original.primalColumnSolution();
+  int numberSets = gubMatrix->numberSets();
+  const double * cost = original.objective();
+  int lastOdd = gubMatrix->firstAvailable();
+  //assert (numberTotalColumns==numberColumns+numberSlacks);
+  int numberRows = original.numberRows();
+  //int numberStaticRows = gubMatrix->numberStaticRows();
+  const int * startSet = gubMatrix->startSets();
+  unsigned char * status = original.statusArray();
+  unsigned char * rowStatus = status+numberColumns;
+  //assert (firstOdd==numberNormal);
+  for (int i=0;i<numberSets;i++) {
+    int iRow = whichRows[i+numberNonGub];
+    original.setRowStatus(iRow,ClpSimplex::atLowerBound);
+  }
+  const int * id = gubMatrix->id();
+  const double * columnLower = gubMatrix->columnLower();
+  const double * columnUpper = gubMatrix->columnUpper();
+  for (int i = 0; i < numberGubColumns; i++) {
+    int iOrig = whichColumns[i+numberNormal];
+    if (iOrig<numberColumns) {
+      if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
+        originalSolution[iOrig] = columnUpper[i];
+        status[iOrig] = 2;
+      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound && columnLower) {
+        originalSolution[iOrig] = columnLower[i];
+        status[iOrig] = 3;
+      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
+        int iSet = gubMatrix->whichSet(i);
+        originalSolution[iOrig] = gubMatrix->keyValue(iSet);
+        status[iOrig] = 1;
+      } else {
+        originalSolution[iOrig] = 0.0;
+        status[iOrig] = 4;
+      }
+    } else {
+      // slack
+      int iSet = iOrig - numberColumns;
+      int iRow = whichRows[iSet+numberNonGub];
+      if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atUpperBound) {
+        original.setRowStatus(iRow,ClpSimplex::atLowerBound);
+      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::atLowerBound) {
+        original.setRowStatus(iRow,ClpSimplex::atUpperBound);
+      } else if (gubMatrix->getDynamicStatus(i) == ClpDynamicMatrix::soloKey) {
+        original.setRowStatus(iRow,ClpSimplex::basic);
+      }
+    }
+  }
+  for (int i = 0; i < numberNormal; i++) {
+    int iOrig = whichColumns[i];
+    ClpSimplex::Status thisStatus = getStatus(i);
+    if (thisStatus == ClpSimplex::basic)
+      status[iOrig] = 1;
+    else if (thisStatus == ClpSimplex::atLowerBound)
+      status[iOrig] = 3;
+    else if (thisStatus == ClpSimplex::atUpperBound)
+      status[iOrig] = 2;
+    else if (thisStatus == ClpSimplex::isFixed)
+      status[iOrig] = 5;
+    else
+      abort();
+    originalSolution[iOrig] = solution[i];
+  }
+  for (int i = numberNormal; i < lastOdd; i++) {
+    int iOrig = whichColumns[id[i-numberNormal] + numberNormal];
+    if (iOrig<numberColumns) {
+      ClpSimplex::Status thisStatus = getStatus(i);
+      if (thisStatus == ClpSimplex::basic)
+        status[iOrig] = 1;
+      else if (thisStatus == ClpSimplex::atLowerBound)
+        status[iOrig] = 3;
+      else if (thisStatus == ClpSimplex::atUpperBound)
+        status[iOrig] = 2;
+      else if (thisStatus == ClpSimplex::isFixed)
+        status[iOrig] = 5;
+      else
+        abort();
+      originalSolution[iOrig] = solution[i];
+    } else {
+      // slack (basic probably)
+      int iSet = iOrig - numberColumns;
+      int iRow = whichRows[iSet+numberNonGub];
+      ClpSimplex::Status thisStatus = getStatus(i);
+      if (thisStatus == ClpSimplex::atLowerBound)
+        thisStatus = ClpSimplex::atUpperBound;
+      else if (thisStatus == ClpSimplex::atUpperBound)
+        thisStatus = ClpSimplex::atLowerBound;
+      original.setRowStatus(iRow,thisStatus);
+    }
+  }
+  for (int i = 0; i < numberNonGub; i++) {
+    int iOrig = whichRows[i];
+    ClpSimplex::Status thisStatus = getRowStatus(i);
+    if (thisStatus == ClpSimplex::basic)
+      rowStatus[iOrig] = 1;
+    else if (thisStatus == ClpSimplex::atLowerBound)
+      rowStatus[iOrig] = 3;
+    else if (thisStatus == ClpSimplex::atUpperBound)
+      rowStatus[iOrig] = 2;
+    else if (thisStatus == ClpSimplex::isFixed)
+      rowStatus[iOrig] = 5;
+    else
+      abort();
+  }
+  int * numberKey = new int [numberRows];
+  memset(numberKey,0,numberRows*sizeof(int));
+  for (int i=0;i<numberSets;i++) {
+    int iRow = whichRows[i+numberNonGub];
+    for (int j=startSet[i];j<startSet[i+1];j++) {
+      int iOrig = whichColumns[j+numberNormal];
+      if (iOrig<numberColumns) {
+        if (original.getColumnStatus(iOrig)==ClpSimplex::basic) {
+          numberKey[iRow]++;
+        }
+      } else {
+        // slack
+        if (original.getRowStatus(iRow)==ClpSimplex::basic) 
+          numberKey[iRow]++;
+      }
+    }
+  }
+  for (int i=0;i<numberSets;i++) {
+    int iRow = whichRows[i+numberNonGub];
+    if (!numberKey[iRow]) {
+      original.setRowStatus(iRow,ClpSimplex::basic);
+    }
+  }
+  delete [] numberKey;
+  double objValue = 0.0;
+  for (int i = 0; i < numberColumns; i++)
+    objValue += cost[i] * originalSolution[i];
+  //printf("objective value is %g\n", objValue);
+}

trunk/Clp/src/ClpSimplexOther.hpp

@@ -190 +190 @@
                       const int * whichRows, const int * whichColumns,
                       int nBound);
+     /** Returns gub version of model or NULL
+         whichRows has to be numberRows
+         whichColumns has to be numberRows+numberColumns */
+     ClpSimplex * gubVersion(int * whichRows, int * whichColumns,
+                             int neededGub,
+                             int factorizationFrequency=50);
+     /// Sets basis from original
+     void setGubBasis(const ClpSimplex &original,const int * whichRows,
+                      const int * whichColumns);
+     /// Restores basis to original
+     void getGubBasis(ClpSimplex &original,const int * whichRows,
+                      const int * whichColumns) const;
      /// Quick try at cleaning up duals if postsolve gets wrong
      void cleanupAfterPostsolve();

trunk/Clp/src/ClpSimplexPrimal.cpp

@@ -1092 +1092 @@
           objectiveWork_ = cost_;
           rowObjectiveWork_ = cost_ + numberColumns_;
-          if (n)
+          if (n||matrix_->type()>=15)
                gutsOfSolution(NULL, NULL, ifValuesPass != 0);
      }
@@ -1118 +1118 @@
                     && firstFree_ < 0) {
                if (handler_->logLevel() == 63)
-                    printf("lastobj %g this %g force %d ", lastObj, thisObj, forceFactorization_);
+                    printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
                int maxFactor = factorization_->maximumPivots();
                if (maxFactor > 10) {
@@ -1130 +1130 @@
                      && firstFree_ < 0) {
                if (handler_->logLevel() == 63)
-                    printf("lastobj3 %g this3 %g `force %d ", lastObj3, thisObj, forceFactorization_);
+                    printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
                int maxFactor = factorization_->maximumPivots();
                if (maxFactor > 10) {

trunk/Clp/src/ClpSolve.cpp

@@ -3025 +3025 @@
      if (matched == (1 << (CLP_PROGRESS - 1)))
           numberMatched = 0;
-     if (numberMatched) {
+     if (numberMatched && model_->clpMatrix()->type() < 15) {
           model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
                     << numberMatched