Changeset 190 for branches/pre

Timestamp:

Aug 5, 2003 4:57:40 AM (18 years ago)

Author:

forrest

Message:

Stuff

Location:

branches/pre

Files:

• ClpMessage.cpp (modified) (1 diff)
• ClpNonLinearCost.cpp (modified) (1 diff)
• ClpPrimalQuadraticDantzig.cpp (modified) (2 diffs)
• ClpSimplex.cpp (modified) (1 diff)
• ClpSimplexPrimal.cpp (modified) (4 diffs)
• ClpSimplexPrimalQuadratic.cpp (modified) (54 diffs)
• include/ClpMessage.hpp (modified) (1 diff)

branches/pre/ClpMessage.cpp

@@ -61 +61 @@
   {CLP_CRASH,28,1,"Crash put %d variables in basis, %d dual infeasibilities"},
   {CLP_END_VALUES_PASS,29,1,"End of values pass after %d iterations"},
+  {CLP_QUADRATIC_BOTH,108,32,"%s %d (%g) and %d (%g) both basic"},
+  {CLP_QUADRATIC_PRIMAL_DETAILS,109,32,"coeff %g, %g, %g - dj %g - deriv zero at %g, sj at %g"},
   {CLP_DUMMY_END,999999,0,""}
 };

branches/pre/ClpNonLinearCost.cpp

@@ -661 +661 @@
   int end = start_[iPivot+1]-1;
   if (!bothWays_) {
-    for (iRange=start; iRange<end;iRange++) {
-      if (value<lower_[iRange+1]+primalTolerance) {
-        // put in better range
-        if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
-          iRange++;
-        break;
+    // If fixed try and get feasible
+    if (lower_[start+1]==lower_[start+2]&&fabs(value-lower_[start+1])<1.001*primalTolerance) {
+      iRange =start+1;
+    } else {
+      for (iRange=start; iRange<end;iRange++) {
+        if (value<=lower_[iRange+1]+primalTolerance) {
+          // put in better range
+          if (value>=lower_[iRange+1]-primalTolerance&&infeasible(iRange)&&iRange==start) 
+            iRange++;
+          break;
+        }
       }
     }

branches/pre/ClpPrimalQuadraticDantzig.cpp

@@ -107 +107 @@
   const double * objective = quadraticInfo_->originalObjective();
   for (iSequence=0;iSequence<originalNumberColumns;iSequence++) {
+    if (model_->flagged(iSequence))
+      continue;
     int jSequence = which[iSequence];
     double value;
@@ -163 +165 @@
   for (jSequence=0;jSequence<originalNumberRows;jSequence++) {
     int iSequence  = jSequence + firstSlack;
+    if (model_->flagged(iSequence))
+      continue;
     int iPi = jSequence+originalNumberColumns;
     // for slacks either pi zero or row at bound

branches/pre/ClpSimplex.cpp

@@ -1081 +1081 @@
     }
   }
-#if 0
-  int status=-99;
-  int * rowIsBasic = new int[numberRows_];
-  int * columnIsBasic = new int[numberColumns_];
-  //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
-  while (status<-98) {
-  
+  if (0)  {
+    static int k=0;
+    printf("start basis\n");
     int i;
-    int numberBasic=0;
-    for (i=0;i<numberRows_;i++) {
-      if (getRowStatus(i) == basic) {
-        rowIsBasic[i]=1;
-        numberBasic++;
-      } else {
-        rowIsBasic[i]=-1;
-      }
-    }
-    for (i=0;i<numberColumns_;i++) {
-      if (getColumnStatus(i) == basic) {
-        columnIsBasic[i]=1;
-        numberBasic++;
-      } else {
-        columnIsBasic[i]=-1;
-      }
-    }
-    assert (numberBasic<=numberRows_);
-    while (status==-99) {
-      status =  factorization_->factorizeOld(this,matrix_,
-                                          numberRows_,numberColumns_,
-                                          rowIsBasic, columnIsBasic,
-                                          0.0);
-      if (status==-99) {
-        // get more memory
-        factorization_->areaFactor(2.0*factorization_->areaFactor());
-      }
-    }
-    if (!status) {
-      // See whether to redo pivot order
-      if (factorization_->needToReorder()||
-          progress_->lastIterationNumber(0)<0) {
-        // do pivot information
-        for (i=0;i<numberRows_;i++) {
-          if (getRowStatus(i) == basic) {
-            pivotVariable_[rowIsBasic[i]]=i+numberColumns_;
-          }
-        }
-        for (i=0;i<numberColumns_;i++) {
-          if (getColumnStatus(i) == basic) {
-            pivotVariable_[columnIsBasic[i]]=i;
-          }
-        }
-      }
-    } else {
-      // leave pivotVariable_ in useful form for cleaning basis
-      for (i=0;i<numberRows_;i++) {
-        pivotVariable_[i]=-1;
-      }
-      for (i=0;i<numberRows_;i++) {
-        if (getRowStatus(i) == basic) {
-          int iPivot = rowIsBasic[i];
-          if (iPivot>=0) 
-            pivotVariable_[iPivot]=i+numberColumns_;
-        }
-      }
-      for (i=0;i<numberColumns_;i++) {
-        if (getColumnStatus(i) == basic) {
-          int iPivot = columnIsBasic[i];
-          if (iPivot>=0) 
-            pivotVariable_[iPivot]=i;
-        }
-      }
-    }
-    if (status==-1) {
-      if (!solveType) {
-        //redo basis - first take ALL columns out
-        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
-          if (getColumnStatus(iColumn)==
-              basic) {
-            // take out
-            if (!valuesPass) {
-              double lower=columnLowerWork_[iColumn];
-              double upper=columnUpperWork_[iColumn];
-              double value=columnActivityWork_[iColumn];
-              if (lower>-largeValue_||upper<largeValue_) {
-                if (fabs(value-lower)<fabs(value-upper)) {
-                  setColumnStatus(iColumn,atLowerBound);
-                  columnActivityWork_[iColumn]=lower;
-                } else {
-                  setColumnStatus(iColumn,atUpperBound);
-                  columnActivityWork_[iColumn]=upper;
-                }
-              } else {
-                setColumnStatus(iColumn,isFree);
-              }
-            } else {
-              setColumnStatus(iColumn,superBasic);
-            }
-          }
-        }
-        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
-              setColumnStatus(iSequence,basic);
-            }
-          } else {
-            // put in slack
-            setRowStatus(iRow,basic);
-          }
-        }
-        // signal repeat
-        status=-99;
-        // set fixed if they are
-        for (iRow=0;iRow<numberRows_;iRow++) {
-          if (getRowStatus(iRow)!=basic ) {
-            if (rowLowerWork_[iRow]==rowUpperWork_[iRow]) {
-              rowActivityWork_[iRow]=rowLowerWork_[iRow];
-              setRowStatus(iRow,isFixed);
-            }
-          }
-        }
-        for (iColumn=0;iColumn<numberColumns_;iColumn++) {
-          if (getColumnStatus(iColumn)!=basic ) {
-            if (columnLowerWork_[iColumn]==columnUpperWork_[iColumn]) {
-              columnActivityWork_[iColumn]=columnLowerWork_[iColumn];
-              setColumnStatus(iColumn,isFixed);
-            }
-          }
-        }
-      }
-    }
-  } 
-  delete [] rowIsBasic;
-  delete [] columnIsBasic;
-#else
+    for (i=0;i<numberRows_;i++)
+      printf ("xx %d %d\n",i,pivotVariable_[i]);
+    for (i=0;i<numberRows_+numberColumns_;i++)
+      if (getColumnStatus(i)==basic)
+        printf ("yy %d basic\n",i);
+    if (k>20)
+      exit(0);
+    k++;
+  }
   int status = factorization_->factorize(this, solveType,valuesPass);
-#endif
   if (status) {
     handler_->message(CLP_SIMPLEX_BADFACTOR,messages_)

branches/pre/ClpSimplexPrimal.cpp

@@ -1251 +1251 @@
     double change = theta*alpha;
     value -= change;
-
+    // But make sure one going out is feasible
     if (change>0.0) {
       // going down
-      if (value<=lower(iPivot)+primalTolerance_) {
+      if (value<lower(iPivot)+primalTolerance_) {
+        if (iPivot==sequenceOut_&&value>lower(iPivot)-1.001*primalTolerance_)
+          value=lower(iPivot);
         double difference = nonLinearCost_->setOne(iPivot,value);
         work[iRow] = difference;
@@ -1267 +1269 @@
     } else {
       // going up
-      if (value>=upper(iPivot)-primalTolerance_) {
+      if (value>upper(iPivot)-primalTolerance_) {
+        if (iPivot==sequenceOut_&&value<upper(iPivot)+1.001*primalTolerance_)
+          value=upper(iPivot);
         double difference = nonLinearCost_->setOne(iPivot,value);
         work[iRow] = difference;
@@ -1685 +1689 @@
              valueOut_<=upperValue+primalTolerance_);
       // may not be exactly at bound and bounds may have changed
+      // Make sure outgoing looks feasible
+      double useValue=valueOut_;
       if (valueOut_<=lowerValue+primalTolerance_) {
         directionOut_=1;
+        useValue= lower_[sequenceOut_];
       } else if (valueOut_>=upperValue-primalTolerance_) {
         directionOut_=-1;
+        useValue= upper_[sequenceOut_];
       } else {
         printf("*** variable wandered off bound %g %g %g!\n",
                lowerValue,valueOut_,upperValue);
         if (valueOut_-lowerValue<=upperValue-valueOut_) {
+          useValue= lower_[sequenceOut_];
           valueOut_=lowerValue;
           directionOut_=1;
         } else {
+          useValue= upper_[sequenceOut_];
           valueOut_=upperValue;
           directionOut_=-1;
@@ -1701 +1711 @@
       }
       solution_[sequenceOut_]=valueOut_;
-      nonLinearCost_->setOne(sequenceOut_,valueOut_);
+      nonLinearCost_->setOne(sequenceOut_,useValue);
     }
     // change cost and bounds on incoming if primal

branches/pre/ClpSimplexPrimalQuadratic.cpp

@@ -563 +563 @@
       if (lastGoodIteration_==numberIterations_)
         factorType=3;
-      if (info->currentPhase())
-        info->setCurrentSolution(solution_);
       // may factorize, checks if problem finished
       statusOfProblemInPrimal(lastCleaned,factorType,progress_,info);
-      if (info->currentPhase())
-        memcpy(solution_,info->currentSolution(),
-               (numberRows_+numberColumns_)*sizeof(double));
       
       // Compute objective function from scratch
@@ -610 +605 @@
             if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
                 ||getColumnStatus(iColumn)==superBasic) {
-              if (getColumnStatus(iColumn)!=basic) {
+              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
                 if (phase!=2) {
                   phase=2;
@@ -636 +631 @@
           if ((value>lower+primalTolerance_&&value<upper-primalTolerance_)
             ||getColumnStatus(iColumn)==superBasic) {
-            if (getColumnStatus(iColumn)!=basic) {
+            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
               if (phase!=2) {
                 phase=2;
@@ -659 +654 @@
 
       // exit if victory declared
-      if (!info->currentPhase()&&primalColumnPivot_->pivotColumn(rowArray_[1],
+      if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1],
                                           rowArray_[2],rowArray_[3],
                                           columnArray_[0],
@@ -693 +688 @@
   checkComplementarity (info,rowArray_[0],rowArray_[1]);
   int crucialSj = info->crucialSj();
+  if (info->crucialSj()>=0) {
+    printf("after inv %d Sj value %g\n",crucialSj,solution_[info->crucialSj()]);
+  }
   int returnCode=-1;
   int phase = info->currentPhase();
   double saveObjective = objectiveValue_;
   int numberXColumns = info->numberXColumns();
+  int numberXRows = info->numberXRows();
   const int * which = info->quadraticSequence();
   const double * pi = solution_+numberXColumns;
@@ -706 +705 @@
   int nextSequenceIn=info->sequenceIn();
   int oldSequenceIn=nextSequenceIn;
+  int saveSequenceIn = sequenceIn_;
   // status stays at -1 while iterating, >=0 finished, -2 to invert
   // status -3 to go to top without an invert
@@ -730 +730 @@
     // Only do if one not already chosen
     bool cleanupIteration;
-    if (phase==2) {
-      // values pass
-      if (nextSequenceIn<0) {
+    if (nextSequenceIn<0) {
+      cleanupIteration=false;
+      if (phase==2) {
+        // values pass
         // get next
         int iColumn;
@@ -741 +742 @@
           double upper = upper_[iColumn];
           if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
-            if (getColumnStatus(iColumn)!=basic) {
+            if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
               nextSequenceIn=iColumn;
               break;
@@ -756 +757 @@
             double upper = upper_[iColumn];
             if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
-              if (getColumnStatus(iColumn)!=basic) {
+              if (getColumnStatus(iColumn)!=basic&&!flagged(iColumn)) {
                 nextSequenceIn=iColumn;
                 break;
@@ -763 +764 @@
           }
         }
-      }
-      oldSequenceIn=nextSequenceIn;
-      sequenceIn_ = nextSequenceIn;
-      cleanupIteration=false;
-      lowerIn_=lower_[sequenceIn_];
-      upperIn_=upper_[sequenceIn_];
-      if (sequenceIn_<numberXColumns) {
-        int jSequence = which[sequenceIn_];
-        if (jSequence>=0) {
-          dualIn_ = solution_[jSequence+numberRows_]; 
-        } else {
-          // need coding
-          dualIn_=cost_[sequenceIn_];
-          int j;
-          for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
-            int iRow = row[j];
-            dualIn_ -= element[j]*pi[iRow];
-          }
-        }
+        oldSequenceIn=nextSequenceIn;
+        saveSequenceIn=sequenceIn_;
+        sequenceIn_ = nextSequenceIn;
       } else {
-        dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
-      }
-      valueIn_=solution_[sequenceIn_];
-      if (dualIn_>0.0)
-        directionIn_ = -1;
-      else 
-        directionIn_ = 1;
-    } else {
-      if (nextSequenceIn<0) {
+        saveSequenceIn=sequenceIn_;
         primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
                      columnArray_[0],columnArray_[1]);
+      }
+    } else {
+      saveSequenceIn=sequenceIn_;
+      sequenceIn_ = nextSequenceIn;
+      if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_) 
         cleanupIteration=false;
-      } else {
-        sequenceIn_ = nextSequenceIn;
+      else
         cleanupIteration=true;
-      }
     }
     pivotRow_=-1;
@@ -811 +791 @@
       // do second half of iteration
       while (chosen>=0) {
+        {
+          // Compute objective function from scratch
+          const CoinPackedMatrix * quadratic = 
+            info->originalModel()->quadraticObjective();
+          const int * columnQuadratic = quadratic->getIndices();
+          const int * columnQuadraticStart = quadratic->getVectorStarts();
+          const int * columnQuadraticLength = quadratic->getVectorLengths();
+          const double * quadraticElement = quadratic->getElements();
+          const double * originalCost = info->originalObjective();
+          int iColumn;
+          objectiveValue_=0.0;
+          double infeasCost=0.0;
+          double linearCost=0.0;
+          for (iColumn=0;iColumn<numberXColumns;iColumn++) {
+            double valueI = solution_[iColumn];
+            linearCost += valueI*originalCost[iColumn];
+            double diff =cost_[iColumn]-originalCost[iColumn];
+            // WE are always feasible so look at low not up!
+            if (diff>0.0) {
+              double above = valueI-lower_[iColumn];
+              assert(above>=0.0);
+              infeasCost += diff*above;
+            } else if (diff<0.0) {
+              double below = upper_[iColumn]-valueI;
+              assert(below>=0.0);
+              infeasCost -= diff*below;
+            }
+            int j;
+            for (j=columnQuadraticStart[iColumn];
+                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
+              int jColumn = columnQuadratic[j];
+              double valueJ = solution_[jColumn];
+              double elementValue = 0.5*quadraticElement[j];
+              objectiveValue_ += valueI*valueJ*elementValue;
+            }
+          }
+          int iRow;
+          for (iRow=0;iRow<numberXRows;iRow++) {
+            int iColumn = iRow + numberColumns_;
+            double diff =cost_[iColumn];
+            double valueI = solution_[iColumn];
+            if (diff>0.0) {
+              double above = valueI-lower_[iColumn];
+              assert(above>=0.0);
+              infeasCost += diff*above;
+            } else if (diff<0.0) {
+              double below = upper_[iColumn]-valueI;
+              assert(below>=0.0);
+              infeasCost -= diff*below;
+            }
+          }
+          objectiveValue_ += linearCost;
+          assert (infeasCost>=0.0);
+          printf("True objective is %g, infeas cost %g, sum %g\n",
+                 objectiveValue_,infeasCost,objectiveValue_+infeasCost);
+        }
         objectiveValue_=saveObjective;
         returnCode=-1;
@@ -841 +877 @@
                 iSjRow = iRow;
                 double d2 = solution_[crucialSj]/tj;
+                if (fabs(solution_[crucialSj])<1.0e-8)
+                  printf("zero crucialSj - pivot out - get right way\n");
+                //if (sequenceIn_>=numberColumns_)
+                //d2=-d2; // Try both ways? - makes no difference!!!
+                //if (sequenceIn_>=numberXColumns&&sequenceIn_<numberColumns_)
+                //d2=-d2; // Try both ways? - makes no difference!!!
                 // see which way to go
                 if (d2>0)
@@ -865 +907 @@
             directionIn_ = 1;
         } else {
+          // not clean up
+          lowerIn_=lower_[sequenceIn_];
+          upperIn_=upper_[sequenceIn_];
+          if (sequenceIn_<numberXColumns) {
+            int jSequence = which[sequenceIn_];
+            if (jSequence>=0) {
+              dualIn_ = solution_[jSequence+numberRows_]; 
+            } else {
+              // need coding
+              dualIn_=cost_[sequenceIn_];
+              int j;
+              for (j=columnStart[sequenceIn_];j<columnStart[sequenceIn_]+columnLength[sequenceIn_]; j++) {
+                int iRow = row[j];
+                dualIn_ -= element[j]*pi[iRow];
+              }
+            }
+          } else {
+            dualIn_ = solution_[sequenceIn_-numberColumns_+numberXColumns];
+          }
+          valueIn_=solution_[sequenceIn_];
+          if (dualIn_>0.0)
+            directionIn_ = -1;
+          else 
+            directionIn_ = 1;
           if (sequenceIn_<numberColumns_) {
             // Set dj as value of slack
@@ -878 +944 @@
           info->setCrucialSj(crucialSj);
         }
+        double oldSj=1.0e30;
+        if (info->crucialSj()>=0&&cleanupIteration)
+          oldSj= solution_[info->crucialSj()];
         // save reduced cost
         //double saveDj = dualIn_;
@@ -886 +955 @@
                              info,
                              cleanupIteration);
+        if (pivotRow_==-1&&phase==2&&fabs(dualIn_)<1.0e-3) {
+          // try other way
+          dualIn_=-dualIn_;
+          directionIn_=-directionIn_;
+          if (info->crucialSj()>=0)
+            setColumnStatus(info->crucialSj(),basic);
+          rowArray_[3]->checkClear();
+          rowArray_[2]->checkClear();
+          rowArray_[0]->checkClear();
+          result=primalRow(rowArray_[1],rowArray_[3],
+                           rowArray_[2],rowArray_[0],
+                           info,
+                           cleanupIteration);
+        }
         saveObjective = objectiveValue_;
         if (pivotRow_>=0) {
@@ -913 +996 @@
           } else if (updateStatus==2) {
             // major error
+            // Reset sequenceIn_
+            // sequenceIn_=saveSequenceIn;
+            pivotRow_=-1;
             // better to have small tolerance even if slower
             factorization_->zeroTolerance(1.0e-15);
@@ -1020 +1106 @@
                  valueOut_<=upperValue+primalTolerance_);
           // may not be exactly at bound and bounds may have changed
+          // Make sure outgoing looks feasible
+          double useValue=valueOut_;
           if (valueOut_<=lowerValue+primalTolerance_) {
             directionOut_=1;
+            useValue= lower_[sequenceOut_];
           } else if (valueOut_>=upperValue-primalTolerance_) {
             directionOut_=-1;
+            useValue= upper_[sequenceOut_];
           } else {
             printf("*** variable wandered off bound %g %g %g!\n",
                    lowerValue,valueOut_,upperValue);
             if (valueOut_-lowerValue<=upperValue-valueOut_) {
+              useValue= lower_[sequenceOut_];
               valueOut_=lowerValue;
               directionOut_=1;
             } else {
+              useValue= upper_[sequenceOut_];
               valueOut_=upperValue;
               directionOut_=-1;
@@ -1036 +1128 @@
           }
           solution_[sequenceOut_]=valueOut_;
-          nonLinearCost_->setOne(sequenceOut_,valueOut_);
+          nonLinearCost_->setOne(sequenceOut_,useValue);
         }
         // change cost and bounds on incoming if primal
         nonLinearCost_->setOne(sequenceIn_,valueIn_); 
         int whatNext=housekeeping(objectiveChange);
+        if (info->crucialSj()>=0) {
+          assert(fabs(oldSj)>= fabs(solution_[info->crucialSj()]));
+          printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
+        }
+        {
+          // Compute objective function from scratch
+          const CoinPackedMatrix * quadratic = 
+            info->originalModel()->quadraticObjective();
+          const int * columnQuadratic = quadratic->getIndices();
+          const int * columnQuadraticStart = quadratic->getVectorStarts();
+          const int * columnQuadraticLength = quadratic->getVectorLengths();
+          const double * quadraticElement = quadratic->getElements();
+          const double * originalCost = info->originalObjective();
+          int iColumn;
+          objectiveValue_=0.0;
+          double infeasCost=0.0;
+          double linearCost=0.0;
+          for (iColumn=0;iColumn<numberXColumns;iColumn++) {
+            double valueI = solution_[iColumn];
+            linearCost += valueI*originalCost[iColumn];
+            double diff =cost_[iColumn]-originalCost[iColumn];
+            // WE are always feasible so look at low not up!
+            if (diff>0.0) {
+              double above = valueI-lower_[iColumn];
+              assert(above>=0.0);
+              infeasCost += diff*above;
+            } else if (diff<0.0) {
+              double below = upper_[iColumn]-valueI;
+              assert(below>=0.0);
+              infeasCost -= diff*below;
+            }
+            int j;
+            for (j=columnQuadraticStart[iColumn];
+                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
+              int jColumn = columnQuadratic[j];
+              double valueJ = solution_[jColumn];
+              double elementValue = 0.5*quadraticElement[j];
+              objectiveValue_ += valueI*valueJ*elementValue;
+            }
+          }
+          int iRow;
+          for (iRow=0;iRow<numberXRows;iRow++) {
+            int iColumn = iRow + numberColumns_;
+            double diff =cost_[iColumn];
+            double valueI = solution_[iColumn];
+            if (diff>0.0) {
+              double above = valueI-lower_[iColumn];
+              assert(above>=0.0);
+              infeasCost += diff*above;
+            } else if (diff<0.0) {
+              double below = upper_[iColumn]-valueI;
+              assert(below>=0.0);
+              infeasCost -= diff*below;
+            }
+          }
+          objectiveValue_ += linearCost;
+          assert (infeasCost>=0.0);
+          printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
+                 objectiveValue_,infeasCost,objectiveValue_+infeasCost);
+        }
         if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
           // really free but going to zero
           setStatus(sequenceOut_,isFree);
+          if (sequenceOut_==info->crucialSj())
+            info->setCrucialSj(-1);
         }
         checkComplementarity (info,rowArray_[2],rowArray_[3]);
@@ -1060 +1214 @@
         cleanupIteration = true;
         // may not be correct on second time
-        if (result) {
+        if (result&&(returnCode==-1||returnCode==-2)) {
           assert(sequenceOut_<numberXColumns||
                  sequenceOut_>=numberColumns_);
+          int crucialSj=info->crucialSj();
           if (sequenceOut_<numberXColumns) {
             chosen =sequenceOut_ + numberRows_; // sj variable
@@ -1070 +1225 @@
             if (iRow<numberRows_-numberXColumns) {
               int iPi = iRow+numberXColumns;
-              printf("pi for row %d is %g\n",
-                     iRow,solution_[iPi]);
+              //printf("pi for row %d is %g\n",
+              //     iRow,solution_[iPi]);
               chosen=iPi;
             } else {
@@ -1077 +1232 @@
               abort();
             }
+          }
+          // But double check as original incoming might have gone out
+          if (chosen==crucialSj) {
+            chosen=-1;
+            break; // means original has gone out
+          }
+          if (returnCode==-2) {
+            // factorization
+            nextSequenceIn = chosen;
+            break;
           }
         } else {
@@ -1212 +1377 @@
   // sj 
   int iSjRow=-1;
+  // sj for incoming 
+  int iSjRow2=-1,crucialSj2=-1;
+  if (sequenceIn_<numberColumns_) {
+    const int * which = info->quadraticSequence();
+    crucialSj2 = which[sequenceIn_];
+    if (crucialSj2>=0)
+      crucialSj2 += numberRows_; // sj2 which should go to 0.0
+  } else {
+    crucialSj2 = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
+  }
   double tj = 0.0;
   // Change in objective will be theta*coeff1 + theta*theta*coeff2
@@ -1234 +1409 @@
       }
     }
+    if (iPivot==crucialSj2)
+      iSjRow2=iRow;
   }
   // Incoming
@@ -1239 +1416 @@
     index2[number2++]=sequenceIn_;
     rhs[sequenceIn_]=way;
-    printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
+    //printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     coeff1 += way*cost_[sequenceIn_];
   } else {
-    printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
+    //printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
     coeffSlack += way*cost_[sequenceIn_];
   }
   rhsArray->setNumElements(number2);
-  if (numberIterations_>=0||cleanupIteration) {
+  if (numberIterations_!=-701) {  //if (numberIterations_>=0||cleanupIteration) {
     for (iIndex=0;iIndex<number2;iIndex++) {
       int iColumn=index2[iIndex];
@@ -1307 +1484 @@
   }
   coeff2 *= 0.5;
-  printf("coefficients %g %g %g - dualIn %g\n",coeff1,coeff2,coeffSlack,dualIn_);
   if (coeffSlack)
     printf("trouble\n");
@@ -1314 +1490 @@
   int accuracyFlag=0;
   if (!cleanupIteration) {
-    assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
-    assert (way*coeff1*dualIn_>=0.0);
+    //assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
+    //assert (way*coeff1*dualIn_>=0.0);
     if (way*coeff1*dualIn_<0.0) {
-      // bad
-      accuracyFlag=2;
-    } else if (fabs(way*coeff1-dualIn_)>1.0e-2*(1.0e-3+fabs(dualIn_))) {
       // bad
       accuracyFlag=2;
@@ -1326 +1499 @@
       accuracyFlag=1;
     }
-    dualIn_ = way*coeff1;
     if (crucialSj>=0) {
-      solution_[crucialSj]=dualIn_;
+      solution_[crucialSj]=way*coeff1;
     }
   } else if (dualIn_<0.0&&way*coeff1>1.0e-2) {
@@ -1345 +1517 @@
   else
     d1 = 0.0;
-  if (fabs(tj)<1.0e-7) {
-    if (sequenceIn_<numberXColumns)
-      printf("column %d is basically linear\n",sequenceIn_);
+  if (fabs(tj)<1.0e-7||!cleanupIteration) {
+    //if (sequenceIn_<numberXColumns)
+      //printf("column %d is basically linear\n",sequenceIn_);
     //assert(!columnQuadraticLength[sequenceIn_]);
   } else {
@@ -1356 +1528 @@
     }
   }
-  printf("derivative zero at %g, sj zero at %g\n",d1,d2);
   if (d1>1.0e10&&d2>1.0e10) {
     // linear variable entering
     // maybe we should have done dual iteration to force sj to 0.0
-    printf("linear variable\n");
-  }
-  maximumMovement = min(maximumMovement,d2);
-  if (d1>=0.0) {
-    maximumMovement = min(maximumMovement,d1);
-    d2 = min(d1,d2);
-  }
-    
+    //printf("linear variable\n");
+  }
+  handler_->message(CLP_QUADRATIC_PRIMAL_DETAILS,messages_)
+    <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
+    <<CoinMessageEol;
+  if (!cleanupIteration)
+    dualIn_ = way*coeff1;
+  double mind1d2=1.0e50;
+  if (cleanupIteration) {
+    // we are only interested in d1 if smaller than d2
+    mind1d2 = min(maximumMovement,d2);
+    //if (d1>1.0e-8&&d1<0.999*d2)
+    //mind1d2=d1;
+  } else {
+    // There is no d2 - d1 refers to crucialSj
+    if (d1>1.0e-5) {
+      mind1d2 = min(maximumMovement,d1);
+      //assert (d1>=0.0);
+    }
+  }
+  maximumMovement = min(maximumMovement,mind1d2);
   rhsArray->clear();
   double tentativeTheta = maximumMovement;
   double upperTheta = maximumMovement;
-
+  bool throwAway=false;
+  if (numberIterations_==700||numberIterations_==694) {
+    printf("Bad iteration coming up after iteration %d\n",numberIterations_);
+  }
 
   for (iIndex=0;iIndex<number;iIndex++) {
@@ -1458 +1645 @@
       maximumSwapped = max(maximumSwapped,numberSwapped);
 
-      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 99999999;
+      double dualCheck = - 2.0*coeff2*tentativeTheta - coeff1 - 0.1*infeasibilityCost_;
       // but make a bit more pessimistic
       dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
@@ -1475 +1662 @@
         memset(spare+numberRemaining,0,
                (saveNumber-numberRemaining)*sizeof(double));
+        double lastThru = totalThru;
         int iTry;
 #define MAXTRY 100
@@ -1539 +1727 @@
             break;
           } else {
-            dualCheck = - 2.0*coeff2*theta_ - coeff1-9999999;
+            dualCheck = - 2.0*coeff2*theta_ - coeff1;
             //dualCheck =0.0;
-            if (totalThru>=dualCheck) {
+            if (totalThru>=dualCheck&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
+              if (!cleanupIteration) {
+                // We can pivot out sj
+                if (upperTheta==maximumMovement) {
+                  printf("figures\n");
+                } else {
+                  if (lastThru>=dualCheck)
+                    printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru);
+                }
+              } else {
+                if (pivotRow_<0&&lastPivotRow<0)
+                  throwAway=true;
+              }
+              break; // no point trying another loop
+            } else if (totalThru>=dualCheck||upperTheta==maximumMovement) {
               break; // no point trying another loop
             } else {
@@ -1554 +1756 @@
               if (bestPivot>bestEverPivot)
                 bestEverPivot=bestPivot;
+              lastThru = totalThru;
             }
           }
@@ -1637 +1840 @@
     dualOut_ = reducedCost(sequenceOut_);
   } else {
+    // If no pivot but bad move then throw away
+    if (throwAway&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
+      assert (pivotRow_<0);
+      accuracyFlag=2;
+    }
     double trueMaximumMovement;
     if (way>0.0) 
@@ -1663 +1871 @@
       }
       // we may still have sj to get rid of
-    } else if (fabs(maximumMovement-d2)<dualTolerance_) {
-      // sj going to zero
-      result=0;
+    } else if (fabs(maximumMovement-mind1d2)<dualTolerance_) {
+      // crucialSj local copy i.e. dj going to zero
+      if (maximumMovement==d2) {
+        // sj going to zero
+        result=0;
+      } else {
+        // incoming dj going to zero
+        assert(!cleanupIteration);
+        if (!cleanupIteration)
+          result=0;
+        iSjRow=iSjRow2;
+        crucialSj = pivotVariable_[iSjRow];
+      }
       assert (pivotRow_<0);
       setColumnStatus(crucialSj,isFree);
@@ -1672 +1890 @@
       // translate to sequence
       sequenceOut_ = pivotVariable_[pivotRow_];
+      assert (sequenceOut_==crucialSj);
       valueOut_ = solution(sequenceOut_);
-      theta_ = d2;
+      theta_ = maximumMovement;
       if (way<0.0) 
         theta_ = - theta_;
@@ -1697 +1916 @@
   // Change in objective will be theta*coeff1 + theta*theta*coeff2
   objectiveValue_ += theta_*coeff1+theta_*theta_*coeff2;
-  printf("New objective value %g\n",objectiveValue_);
   {
     int iColumn;
@@ -1722 +1940 @@
       }
     }
-    printf("Objective value %g\n",objectiveValue());
   }
   if (accuracyFlag<2) {
@@ -1737 +1954 @@
                                                    ClpQuadraticInfo * info)
 {
+  if (info->currentPhase())
+    info->setCurrentSolution(solution_);
+  //info->saveStatus();
+  
   if (type==2) {
     // trouble - restore solution
@@ -1761 +1982 @@
     if (type) {
       // is factorization okay?
+      int numberPivots = factorization_->pivots();
       if (internalFactorize(1)) {
         if (solveType_==2) {
@@ -1767 +1989 @@
           return;
         }
-
+        numberIterations_ = progress_->lastIterationNumber(0);
         // no - restore previous basis
         assert (type==1);
@@ -1779 +2001 @@
         assert (internalFactorize(1)==0);
         info->restoreStatus();
+        // flag incoming
+        if (numberPivots==1) {
+          int crucialSj = info->crucialSj();
+          if (crucialSj<0) {
+            setFlagged(sequenceIn_);
+          } else {
+            int iSequence;
+            int numberXColumns = info->numberXColumns();
+            if (crucialSj<numberRows_) {
+              // row
+              iSequence = crucialSj-numberXColumns+numberColumns_;
+            } else {
+              // column
+              iSequence = crucialSj-numberRows_;
+            }
+            if (getColumnStatus(iSequence)!=basic) {
+              setFlagged(iSequence);
+              info->setSequenceIn(-1);
+            } else {
+              printf("**** %d is basic ?\n",iSequence);
+            }
+          }
+        }
         changeMade_++; // say change made
       }
@@ -1789 +2034 @@
   // put back original costs and then check
   createRim(4);
+  if (info->currentPhase()) {
+    memcpy(solution_,info->currentSolution(),
+           (numberRows_+numberColumns_)*sizeof(double));
+  }
   gutsOfSolution(NULL,NULL);
+  if (info->currentPhase()) {
+    memcpy(solution_,info->currentSolution(),
+           (numberRows_+numberColumns_)*sizeof(double));
+    nonLinearCost_->checkInfeasibilities(primalTolerance_);
+  }
   checkComplementarity(info,rowArray_[2],rowArray_[3]);
   // Double check reduced costs if no action
@@ -1799 +2053 @@
   }
   // Check if looping
+  // Take out for now
   int loop;
-  if (type!=2) 
+  if (type!=2&&0) 
     loop = progress->looping();
   else
@@ -1852 +2107 @@
       // put back original costs
       createRim(4);
+      // put all non-basic variables to bounds
+      int iSequence;
+      for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
+        Status status = getStatus(iSequence);
+        if (status==atLowerBound||status==isFixed)
+          solution_[iSequence]=lower_[iSequence];
+        else if (status==atUpperBound)
+          solution_[iSequence]=upper_[iSequence];
+      }
       nonLinearCost_->checkInfeasibilities(primalTolerance_);
       gutsOfSolution(NULL,NULL);
+      nonLinearCost_->checkInfeasibilities(primalTolerance_);
       checkComplementarity(info,rowArray_[2],rowArray_[3]);
 
@@ -1865 +2130 @@
         problemStatus_ = -1;
         infeasibilityCost_=saveWeight;
+        nonLinearCost_->checkInfeasibilities(primalTolerance_);
+        checkComplementarity(info,rowArray_[2],rowArray_[3]);
       } else {
         nonLinearCost_=NULL;
@@ -1888 +2155 @@
           infeasibilityCost_=0.0;
           createRim(4);
+          // put all non-basic variables to bounds
+          int iSequence;
+          for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
+            Status status = getStatus(iSequence);
+            if (status==atLowerBound||status==isFixed)
+              solution_[iSequence]=lower_[iSequence];
+            else if (status==atUpperBound)
+              solution_[iSequence]=upper_[iSequence];
+          }
           nonLinearCost_->checkInfeasibilities(primalTolerance_);
           gutsOfSolution(NULL,NULL);
@@ -1940 +2216 @@
           // put back original costs and then check
           createRim(4);
+          // put all non-basic variables to bounds
+          int iSequence;
+          for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
+            Status status = getStatus(iSequence);
+            if (status==atLowerBound||status==isFixed)
+              solution_[iSequence]=lower_[iSequence];
+            else if (status==atUpperBound)
+              solution_[iSequence]=upper_[iSequence];
+          }
           nonLinearCost_->checkInfeasibilities(oldTolerance);
           gutsOfSolution(NULL,NULL);
@@ -2054 +2339 @@
     // zero out basic values
     int iRow;
-    double * save = new double [numberRows_];
     for (iRow=0;iRow<numberRows_;iRow++) {
       int iPivot=pivotVariable_[iRow];
-      save[iRow]=solution_[iPivot];
+      // save
+      dj_[iPivot]=solution_[iPivot];
       solution_[iPivot] = 0.0;
     }
@@ -2103 +2388 @@
       int iSequence  = iRow + numberColumns_;
       double value = cost_[iSequence];
-      if (value) {
-        assert(getRowStatus(iRow)==basic);
+      // Its possible a fixed vector may have had this set
+      if (value&&getRowStatus(iRow)==basic) {
         // add in pi column
         matrix_->add(this,array1,iRow+numberXColumns,value);
@@ -2116 +2401 @@
       int iRow = index[k];
       double value=modifiedCost[iRow];
-      int iSequence = pivotVariable_[iRow];
-      if (iSequence>=numberXColumns&&iSequence<numberColumns_)
-        solution_[iSequence] = value;
+      int iPivot = pivotVariable_[iRow];
+      if (iPivot>=numberXColumns&&iPivot<numberColumns_)
+        solution_[iPivot] = value;
       else
-        solution_[iSequence] = save[iRow];
-      //solution_[iSequence] = value;
+        solution_[iPivot] = dj_[iPivot];
+      //solution_[iPivot] = value;
       
     }
-    delete [] save;
     array1->clear();
     // Now modify pi values by slack costs to make djs
@@ -2144 +2428 @@
     for (k=0;k<numberRows_;k++) {
       if (k<numberXRows) {
-        if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
-          printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
+        //if (fabs(modifiedCost[k]-rowSolution[k])>1.0e-3*(1.0+fabs(rowSolution[k])))
+        //printf("bad row %d %g %g\n",k,modifiedCost[k],rowSolution[k]);
       } else {
         rowSolution[k]=modifiedCost[k];
@@ -2204 +2488 @@
         jSequence += start;
         if (getColumnStatus(jSequence)==basic)
-          printf("Struct %d (%g) and %d (%g) both basic\n",
-                 iSequence,solution_[iSequence],jSequence,solution_[jSequence]);
+          handler_->message(CLP_QUADRATIC_BOTH,messages_)
+            <<"Structural"<<iSequence
+            <<solution_[iSequence]<<jSequence<<solution_[jSequence]
+            <<CoinMessageEol;
       }
       break;
@@ -2850 +3136 @@
   int numberRows = numberXRows_+numberQuadraticColumns_;
   int size = numberRows+numberColumns;
-  if (!validSolution_) 
-    validSolution_ = new double [size];
-  memcpy(validSolution_,currentSolution_,size*sizeof(double));
+  if (currentSolution_) {
+    if (!validSolution_) 
+      validSolution_ = new double [size];
+    memcpy(validSolution_,currentSolution_,size*sizeof(double));
+  }
 }
 // Restore previous

branches/pre/include/ClpMessage.hpp

@@ -61 +61 @@
   CLP_CRASH,
   CLP_END_VALUES_PASS,
+  CLP_QUADRATIC_BOTH,
+  CLP_QUADRATIC_PRIMAL_DETAILS,
   CLP_DUMMY_END
 };