Changeset 1677 for trunk

Timestamp:

Jan 21, 2011 1:01:03 PM (9 years ago)

Author:

forrest

Message:

gub

Location:

trunk/Clp/src

Files:

• ClpDynamicMatrix.cpp (modified) (11 diffs)
• ClpMessage.cpp (modified) (1 diff)
• ClpMessage.hpp (modified) (1 diff)
• ClpNonLinearCost.cpp (modified) (1 diff)
• ClpSimplexOther.cpp (modified) (7 diffs)
• OsiClp/OsiClpSolverInterface.cpp (modified) (4 diffs)
• OsiClp/OsiClpSolverInterface.hpp (modified) (2 diffs)

trunk/Clp/src/ClpDynamicMatrix.cpp

@@ -94 +94 @@
      lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
      upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
-     status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
+     status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_+4*sizeof(int));
      model_ = rhs.model_;
      sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
@@ -237 +237 @@
      flags_ &= ~1;
      // resize model (matrix stays same)
-     int newRowSize = numberRows + CoinMin(numberSets_, CoinMax(frequency, numberRows)) + 1;
+     int newRowSize = numberRows + CoinMin(numberSets_, frequency+numberRows) + 1;
      model->resize(newRowSize, numberNeeded);
      for (i = numberRows; i < newRowSize; i++)
@@ -259 +259 @@
      keyVariable_ = new int[numberSets_];
      if (status) {
-          status_ = ClpCopyOfArray(status, 2*numberSets_);
+          status_ = ClpCopyOfArray(status, 2*numberSets_+4*sizeof(int));
           assert (dynamicStatus);
           dynamicStatus_ = ClpCopyOfArray(dynamicStatus, 2*numberGubColumns_);
      } else {
-          status_ = new unsigned char [2*numberSets_];
+          status_ = new unsigned char [2*numberSets_+4*sizeof(int)];
           memset(status_, 0, numberSets_);
           int i;
@@ -356 +356 @@
           lowerSet_ = ClpCopyOfArray(rhs.lowerSet_, numberSets_);
           upperSet_ = ClpCopyOfArray(rhs.upperSet_, numberSets_);
-          status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_);
+          status_ = ClpCopyOfArray(rhs.status_, 2*numberSets_+4*sizeof(int));
           model_ = rhs.model_;
           sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
@@ -865 +865 @@
           if (lowerSet_[iSet] == upperSet_[iSet])
                setStatus(iSet, ClpSimplex::isFixed);
+#if 0
           if (getStatus(iSet) != model->getStatus(sequenceOut))
                printf("** set %d status %d, var status %d\n", iSet,
                       getStatus(iSet), model->getStatus(sequenceOut));
+#endif
      }
      ClpMatrixBase::updatePivot(model, oldInValue, oldOutValue);
@@ -1044 +1046 @@
 {
      int returnCode = 0;
+#if 0 //ndef NDEBUG
+     {
+       int numberColumns = model->numberColumns();
+       int numberRows = model->numberRows();
+       int * pivotVariable = model->pivotVariable();
+       if (pivotVariable&&model->numberIterations()) {
+         for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
+           assert (pivotVariable[i]==i+numberColumns);
+         }
+       }
+     }
+#endif
      switch (mode) {
           // Fill in pivotVariable
@@ -1087 +1101 @@
      case 5: {
        memcpy(status_+numberSets_,status_,numberSets_);
+       memcpy(status_+2*numberSets_,&numberActiveSets_,sizeof(int));
        memcpy(dynamicStatus_+maximumGubColumns_,
               dynamicStatus_,maximumGubColumns_);
@@ -1094 +1109 @@
      case 6: {
        memcpy(status_,status_+numberSets_,numberSets_);
+       memcpy(&numberActiveSets_,status_+2*numberSets_,sizeof(int));
        memcpy(dynamicStatus_,dynamicStatus_+maximumGubColumns_,
               maximumGubColumns_);
@@ -1266 +1282 @@
           return 1;
      }
+#ifndef NDEBUG
+     {
+       int numberColumns = model->numberColumns();
+       int numberRows = model->numberRows();
+       int * pivotVariable = model->pivotVariable();
+       for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
+         assert (pivotVariable[i]==i+numberColumns);
+       }
+     }
+#endif
      // lookup array
      int * active = new int [numberActiveSets_];
@@ -1479 +1505 @@
                }
           }
+          if (iPut<numberStaticRows_+numberActiveSets_) {
+            printf("lost %d sets\n",
+                   numberStaticRows_+numberActiveSets_-iPut);
+            iPut = numberStaticRows_+numberActiveSets_;
+          }
           for (i = numberActiveSets_; i < currentNumberActiveSets; i++) {
                pivotVariable[iPut++] = i + base3;
           }
-          assert (iPut == numberRows);
+          //assert (iPut == numberRows);
      }
 #ifdef CLP_DEBUG
@@ -2288 +2319 @@
           keyVariable_[iSet] = whichKey;
      }
+     // clean up pivotVariable
+     int numberColumns = model_->numberColumns();
+     int numberRows = model_->numberRows();
+     int * pivotVariable = model_->pivotVariable();
+     if (pivotVariable) {
+       for (int i=0; i<numberStaticRows_+numberActiveSets_;i++) {
+         if (model_->getRowStatus(i)!=ClpSimplex::basic)
+           pivotVariable[i]=-1;
+         else
+           pivotVariable[i]=numberColumns+i;
+       }
+       for (int i=numberStaticRows_+numberActiveSets_;i<numberRows;i++) {
+         pivotVariable[i]=i+numberColumns;
+       }
+       int put=-1;
+       for (int i=0;i<numberColumns;i++) {
+         if (model_->getColumnStatus(i)==ClpSimplex::basic) {
+           while(put<numberRows) {
+             put++;
+             if (pivotVariable[put]==-1) {
+               pivotVariable[put]=i;
+               break;
+             }
+           }
+         }
+       }
+       for (int i=CoinMax(put,0);i<numberRows;i++) {
+         if (pivotVariable[i]==-1) 
+           pivotVariable[i]=i+numberColumns;
+       }
+     }
+     if (rhsOffset_) {
+       double * cost = model_->costRegion();
+       double * columnLower = model_->lowerRegion();
+       double * columnUpper = model_->upperRegion();
+       double * solution = model_->solutionRegion();
+       int numberRows = model_->numberRows();
+       for (int i = numberActiveSets_; i < numberRows-numberStaticRows_; i++) {
+         int iSequence = i + numberStaticRows_ + numberColumns;
+         solution[iSequence] = 0.0;
+         columnLower[iSequence] = -COIN_DBL_MAX;
+         columnUpper[iSequence] = COIN_DBL_MAX;
+         cost[iSequence] = 0.0;
+         model_->nonLinearCost()->setOne(iSequence, solution[iSequence],
+                                        columnLower[iSequence],
+                                        columnUpper[iSequence], 0.0);
+         model_->setStatus(iSequence, ClpSimplex::basic);
+         rhsOffset_[i+numberStaticRows_] = 0.0;
+       }
+#if 0
+       for (int i=0;i<numberStaticRows_;i++)
+         printf("%d offset %g\n",
+                i,rhsOffset_[i]);
+#endif
+     }
      numberActiveColumns_ = firstAvailable_;
+#if 0
+     for (iSet = 0; iSet < numberSets_; iSet++) {
+       for (int j=startSet_[iSet];j<startSet_[iSet+1];j++) {
+         if (getDynamicStatus(j)==soloKey)
+           printf("%d in set %d is solo key\n",j,iSet);
+         else if (getDynamicStatus(j)==inSmall)
+           printf("%d in set %d is in small\n",j,iSet);
+       }
+     }
+#endif
      return;
 }

trunk/Clp/src/ClpMessage.cpp

@@ -104 +104 @@
      {CLP_PARAMETRICS_STATS2, 63, 2, "Theta %g - objective %g, %s in, %s out"},
      {CLP_GENERAL, 1000, 1, "%s"},
+     {CLP_GENERAL2, 1001, 2, "%s"},
      {CLP_DUMMY_END, 999999, 0, ""}
 };

trunk/Clp/src/ClpMessage.hpp

@@ -105 +105 @@
      CLP_PARAMETRICS_STATS2,
      CLP_GENERAL,
+     CLP_GENERAL2,
      CLP_DUMMY_END
 };

trunk/Clp/src/ClpNonLinearCost.cpp

@@ -91 +91 @@
           }
      }
-     if (allZero)
+     if (allZero&&model_->clpMatrix()->type()<15)
           model_->setInfeasibilityCost(1.0);
      double infeasibilityCost = model_->infeasibilityCost();

trunk/Clp/src/ClpSimplexOther.cpp

@@ -4173 +4173 @@
   if (numberGub >= neededGub) { 
     sprintf(message,"%d gub rows", numberGub);
-    handler_->message(CLP_GENERAL, messages_)
+    handler_->message(CLP_GENERAL2, messages_)
       << message << CoinMessageEol;
     int numberNormal = 0;
@@ -4183 +4183 @@
     if (!numberNormal) {
       sprintf(message,"Putting back one gub row to make non-empty");
-      handler_->message(CLP_GENERAL, messages_)
+      handler_->message(CLP_GENERAL2, messages_)
         << message << CoinMessageEol;
       rowIsGub[smallestGubRow]=-1;
@@ -4399 +4399 @@
     sprintf(message,"** Before adding matrix there are %d rows and %d columns",
            model2->numberRows(), model2->numberColumns());
-    handler_->message(CLP_GENERAL, messages_)
+    handler_->message(CLP_GENERAL2, messages_)
       << message << CoinMessageEol;
     delete [] scaleArray;
@@ -4439 +4439 @@
             "** While after adding matrix there are %d rows and %d columns",
             model2->numberRows(), model2->numberColumns());
-    handler_->message(CLP_GENERAL, messages_)
+    handler_->message(CLP_GENERAL2, messages_)
       << message << CoinMessageEol;
     model2->setSpecialOptions(4);    // exactly to bound
@@ -4471 +4471 @@
   double * solution = primalColumnSolution();
   const double * originalSolution = original.primalColumnSolution();
+  const double * upperSet = gubMatrix->upperSet();
   int numberSets = gubMatrix->numberSets();
   const int * startSet = gubMatrix->startSets();
@@ -4477 +4478 @@
   for (int i=0;i<numberSets;i++) {
     for (int j=startSet[i];j<startSet[i+1];j++) {
+      gubMatrix->setDynamicStatus(j,ClpDynamicMatrix::atLowerBound);
       int iColumn = whichColumns[j+numberNormal];
       if (iColumn<numberColumns)
@@ -4543 +4545 @@
     int iRow = whichRows[numberNonGub+i];
     if (!numberKey[iRow]) {
-      if (original.getRowStatus(iRow)==ClpSimplex::basic) {
+      double upper = upperSet[i]-1.0e-7;
+      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;
-          }
+      // 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];
+          if (value>upper)
+            gubMatrix->setStatus(i,ClpSimplex::atLowerBound);
+        } else {
+          // slack - take value as 0.0 as will win on length
+          value=0.0;
         }
-        assert(chosen>=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);
+      if (gubMatrix->getStatus(i)!=ClpSimplex::basic) {
         // set as key
         for (int j=startSet[i];j<startSet[i+1];j++) {

trunk/Clp/src/OsiClp/OsiClpSolverInterface.cpp

@@ -57 +57 @@
   }
 #endif
+  if ((specialOptions_&2097152)!=0) {
+    resolveGub((9*modelPtr_->numberRows())/10);
+    return;
+  }
   bool deleteSolver;
   ClpSimplex * solver;
@@ -783 +787 @@
     return;
   }
+  if ((specialOptions_&2097152)!=0) {
+    resolveGub((9*modelPtr_->numberRows())/10);
+    return;
+  }
   //void pclp(char *);
   //pclp("res");
@@ -1269 +1277 @@
     modelPtr_->whatsChanged_ &= ~0xffff;
   modelPtr_->whatsChanged_ |= 0x30000;
+}
+#include "ClpSimplexOther.hpp"
+// Resolve an LP relaxation after problem modification (try GUB)
+void 
+OsiClpSolverInterface::resolveGub(int needed)
+{
+  bool takeHint;
+  OsiHintStrength strength;
+  // Switch off printing if asked to
+  bool gotHint = (getHintParam(OsiDoReducePrint,takeHint,strength));
+  assert (gotHint);
+  int saveMessageLevel=modelPtr_->logLevel();
+  if (strength!=OsiHintIgnore&&takeHint) {
+    int messageLevel=messageHandler()->logLevel();
+    if (messageLevel>0)
+      modelPtr_->messageHandler()->setLogLevel(messageLevel-1);
+    else
+      modelPtr_->messageHandler()->setLogLevel(0);
+  }
+  setBasis(basis_,modelPtr_);
+  // find gub
+  int numberRows = modelPtr_->numberRows();
+  int * which = new int[numberRows];
+  int numberColumns = modelPtr_->numberColumns();
+  int * whichC = new int[numberColumns+numberRows];
+  ClpSimplex * model2 = 
+    static_cast<ClpSimplexOther *> (modelPtr_)->gubVersion(which,whichC,
+                                                           needed);
+  if (model2) {
+    // move in solution
+    static_cast<ClpSimplexOther *> (model2)->setGubBasis(*modelPtr_,
+                                                         which,whichC);
+    model2->setLogLevel(CoinMin(1,model2->logLevel()));
+    ClpPrimalColumnSteepest steepest(5);
+    model2->setPrimalColumnPivotAlgorithm(steepest);
+    double time1 = CoinCpuTime();
+    model2->primal();
+    //printf("Primal took %g seconds\n",CoinCpuTime()-time1);
+    static_cast<ClpSimplexOther *> (model2)->getGubBasis(*modelPtr_,
+                                                         which,whichC);
+    int totalIterations = model2->numberIterations();
+    //modelPtr_->setLogLevel(63);
+    modelPtr_->primal(1);
+    modelPtr_->setNumberIterations(totalIterations+modelPtr_->numberIterations());
+  } else {
+    modelPtr_->dual();
+  }
+  basis_ = getBasis(modelPtr_);
+  modelPtr_->messageHandler()->setLogLevel(saveMessageLevel);
 }
 // Sort of lexicographic resolve
@@ -6276 +6333 @@
       }
       // set normal
-      specialOptions_ &= (2047+3*8192+15*65536);
+      specialOptions_ &= (2047|3*8192|15*65536|2097152);
       if (otherInformation!=NULL) {
         int * array = static_cast<int *> (otherInformation);

trunk/Clp/src/OsiClp/OsiClpSolverInterface.hpp

@@ -52 +52 @@
   /// Resolve an LP relaxation after problem modification
   virtual void resolve();
+  
+  /// Resolve an LP relaxation after problem modification (try GUB)
+  virtual void resolveGub(int needed);
   
   /// Invoke solver's built-in enumeration algorithm
@@ -1369 +1372 @@
       524288 Fake objective and 0-1
       1048576 Don't recompute ray after crunch
+      2097152 
   */
   mutable unsigned int specialOptions_;