Changeset 196 for branches/pre

Timestamp:

Aug 13, 2003 6:43:02 AM (18 years ago)

Author:

forrest

Message:

64 bit, Sbb etc

Location:

branches/pre

Files:

• ClpDualRowSteepest.cpp (modified) (4 diffs)
• ClpFactorization.cpp (modified) (1 diff)
• ClpModel.cpp (modified) (12 diffs)
• ClpPrimalColumnDantzig.cpp (modified) (1 diff)
• ClpPrimalColumnSteepest.cpp (modified) (5 diffs)
• ClpPrimalQuadraticDantzig.cpp (modified) (6 diffs)
• ClpSimplex.cpp (modified) (3 diffs)
• ClpSimplexDual.cpp (modified) (17 diffs)
• ClpSimplexPrimalQuadratic.cpp (modified) (61 diffs)
• Makefile.Clp (modified) (2 diffs)
• Test/ClpMain.cpp (modified) (2 diffs)
• Test/unitTest.cpp (modified) (1 diff)
• include/ClpModel.hpp (modified) (2 diffs)
• include/ClpSimplex.hpp (modified) (1 diff)
• include/ClpSimplexDual.hpp (modified) (1 diff)
• include/ClpSimplexPrimalQuadratic.hpp (modified) (7 diffs)

branches/pre/ClpDualRowSteepest.cpp

@@ -263 +263 @@
     double * array = alternateWeights_->denseVector();
     int * which = alternateWeights_->getIndices();
+    int i;
     for (i=0;i<numberRows;i++) {
       double value=0.0;
@@ -280 +281 @@
       if (fabs(weights_[i]-value)>1.0e-4)
         printf("%d old %g, true %g\n",i,weights_[i],value);
+      //else 
+      //printf("%d matches %g\n",i,value);
     }
     delete temp;
@@ -312 +315 @@
     }
     spare->setNumElements ( numberNonZero );
+    // Only one array active as already permuted
     model_->factorization()->updateColumn(spare,spare,true);
     // permute back
@@ -317 +321 @@
     // alternateWeights_ should still be empty
     int pivotRow = model_->pivotRow();
+#ifdef CLP_DEBUG
+    if ( model_->logLevel (  ) >4  && 
+         fabs(norm-weights_[pivotRow])>1.0e-3*(1.0+norm)) 
+      printf("on row %d, true weight %g, old %g\n",
+             pivotRow,sqrt(norm),sqrt(weights_[pivotRow]));
+#endif
     // could re-initialize here (could be expensive)
     norm /= model_->alpha() * model_->alpha();

branches/pre/ClpFactorization.cpp

@@ -450 +450 @@
 {
 #ifdef CLP_DEBUG
-  regionSparse->checkClear();
+  if (!noPermute)
+    regionSparse->checkClear();
 #endif
   if (!networkBasis_) {

branches/pre/ClpModel.cpp

@@ -20 +20 @@
 #include "ClpMessage.hpp"
 #include "ClpLinearObjective.hpp"
+#include "ClpQuadraticObjective.hpp"
 
 //#############################################################################
@@ -42 +43 @@
   matrix_(NULL),
   rowCopy_(NULL),
-  quadraticObjective_(NULL),
   ray_(NULL),
   status_(NULL),
@@ -106 +106 @@
   delete rowCopy_;
   rowCopy_=NULL;
-  delete quadraticObjective_;
-  quadraticObjective_ = NULL;
   delete [] ray_;
   ray_ = NULL;
@@ -373 +371 @@
       rowCopy_=NULL;
     }
-    if (rhs.quadraticObjective_) {
-      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
-    } else {
-      quadraticObjective_=NULL;
-    }
     matrix_=NULL;
     if (rhs.matrix_) {
@@ -395 +388 @@
     matrix_ = rhs.matrix_;
     rowCopy_ = NULL;
-    quadraticObjective_ = rhs.quadraticObjective_;
     ray_ = rhs.ray_;
     lengthNames_ = 0;
@@ -441 +433 @@
   matrix_ = NULL;
   rowCopy_ = NULL;
-  quadraticObjective_=NULL,
   delete [] otherModel.ray_;
   otherModel.ray_ = ray_;
@@ -677 +668 @@
       which[i-newNumberColumns]=i;
     matrix_->deleteCols(numberColumns_-newNumberColumns,which);
-    if (quadraticObjective_) {
-      quadraticObjective_->deleteCols(numberColumns_-newNumberColumns,which);
-      quadraticObjective_->deleteRows(numberColumns_-newNumberColumns,which);
-    }
     delete [] which;
   }
@@ -752 +739 @@
   matrix_->deleteCols(number,which);
   //matrix_->removeGaps();
-  if (quadraticObjective_) {
-    quadraticObjective_->deleteCols(number,which);
-    quadraticObjective_->deleteRows(number,which);
-  }
   // status
   if (status_) {
@@ -1209 +1192 @@
 {
   assert (numberColumns==numberColumns_);
-  quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
-                                             start[numberColumns],element,column,start,NULL);
+  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
+  double offset;
+  ClpObjective * obj = new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns,
+                                             start,column,element);
+  delete objective_;
+  objective_ = obj;
+
 }
 void 
@@ -1216 +1204 @@
 {
   assert (matrix.getNumCols()==numberColumns_);
-  quadraticObjective_ = new CoinPackedMatrix(matrix);
+  assert ((dynamic_cast< ClpLinearObjective*>(objective_)));
+  double offset;
+  ClpQuadraticObjective * obj = 
+    new ClpQuadraticObjective(objective_->gradient(NULL,offset),numberColumns_,
+                                                 NULL,NULL,NULL);
+  delete objective_;
+  objective_ = obj;
+  obj->loadQuadraticObjective(matrix);
 }
 // Get rid of quadratic objective
@@ -1222 +1217 @@
 ClpModel::deleteQuadraticObjective()
 {
-  delete quadraticObjective_;
-  quadraticObjective_ = NULL;
+  ClpQuadraticObjective * obj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
+  if (obj)
+    obj->deleteQuadraticObjective();
 }
 void 
@@ -1417 +1413 @@
     rowCopy_=NULL;
   }
-  if (rhs->quadraticObjective_) {
-    quadraticObjective_ = new CoinPackedMatrix(*rhs->quadraticObjective_,
-                                               numberColumns,whichColumn,
-                                               numberColumns,whichColumn);
-  } else {
-    quadraticObjective_=NULL;
-  }
   matrix_=NULL;
   if (rhs->matrix_) {

branches/pre/ClpPrimalColumnDantzig.cpp

@@ -85 +85 @@
     anyUpdates=false;
   }
-
+  if (!updates->getNumElements())
+    anyUpdates=false; // in case dj 0.0 (for quadratic)
   if (anyUpdates) {
     model_->factorization()->updateColumnTranspose(spareRow2,updates);

branches/pre/ClpPrimalColumnSteepest.cpp

@@ -414 +414 @@
   }
 
-#ifdef RECOMPUTE_ALL_DJS
-  for (iSection=0;iSection<2;iSection++) {
+  // If in quadratic re-compute all
+  if (model_->algorithm()==2) {
+    for (iSection=0;iSection<2;iSection++) {
       
-    reducedCost=model_->djRegion(iSection);
-    int addSequence;
-    int iSequence;
+      reducedCost=model_->djRegion(iSection);
+      int addSequence;
+      int iSequence;
       
-    if (!iSection) {
-      number = model_->numberRows();
-      addSequence = model_->numberColumns();
-    } else {
-      number = model_->numberColumns();
-      addSequence = 0;
-    }
-
-    if (!model_->nonLinearCost()->lookBothWays()) {
-      for (iSequence=0;iSequence<number;iSequence++) {
-        double value = reducedCost[iSequence];
-        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
-
-        switch(status) {
+      if (!iSection) {
+        number = model_->numberRows();
+        addSequence = model_->numberColumns();
+      } else {
+        number = model_->numberColumns();
+        addSequence = 0;
+      }
+      
+      if (!model_->nonLinearCost()->lookBothWays()) {
+        for (iSequence=0;iSequence<number;iSequence++) {
+          double value = reducedCost[iSequence];
+          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
           
-        case ClpSimplex::basic:
-          infeasible_->zero(iSequence+addSequence);
-        case ClpSimplex::isFixed:
-          break;
-        case ClpSimplex::isFree:
-          if (fabs(value)>tolerance) {
-            // we are going to bias towards free (but only if reasonable)
-            value *= FREE_BIAS;
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
+          switch(status) {
+            
+          case ClpSimplex::basic:
             infeasible_->zero(iSequence+addSequence);
-          }
-          break;
-        case ClpSimplex::atUpperBound:
-          if (value>tolerance) {
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
-            infeasible_->zero(iSequence+addSequence);
-          }
-          break;
-        case ClpSimplex::atLowerBound:
-          if (value<-tolerance) {
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
-            infeasible_->zero(iSequence+addSequence);
-          }
-        }
-      }
-    } else {
-      // we can go both ways
-      ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
-      for (iSequence=0;iSequence<number;iSequence++) {
-        double value = reducedCost[iSequence];
-        ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
-
-        switch(status) {
-          
-        case ClpSimplex::basic:
-          infeasible_->zero(iSequence+addSequence);
-        case ClpSimplex::isFixed:
-          break;
-        case ClpSimplex::isFree:
-          if (fabs(value)>tolerance) {
-            // we are going to bias towards free (but only if reasonable)
-            value *= FREE_BIAS;
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
-            infeasible_->zero(iSequence+addSequence);
-          }
-          break;
-        case ClpSimplex::atUpperBound:
-          if (value>tolerance) {
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
-            // look other way - change up should be negative
-            value -= nonLinear->changeUpInCost(iSequence+addSequence);
-            if (value>-tolerance) {
-              infeasible_->zero(iSequence+addSequence);
-            } else {
+          case ClpSimplex::isFixed:
+            break;
+          case ClpSimplex::isFree:
+          case ClpSimplex::superBasic:
+            if (fabs(value)>tolerance) {
+              // we are going to bias towards free (but only if reasonable)
+              value *= FREE_BIAS;
               // store square in list
               if (infeas[iSequence+addSequence])
@@ -520 +451 @@
               else
                 infeasible_->quickAdd(iSequence+addSequence,value*value);
-            }
-          }
-          break;
-        case ClpSimplex::atLowerBound:
-          if (value<-tolerance) {
-            // store square in list
-            if (infeas[iSequence+addSequence])
-              infeas[iSequence+addSequence] = value*value; // already there
-            else
-              infeasible_->quickAdd(iSequence+addSequence,value*value);
-          } else {
-            // look other way - change down should be positive
-            value -= nonLinear->changeDownInCost(iSequence+addSequence);
-            if (value<tolerance) {
+            } else {
               infeasible_->zero(iSequence+addSequence);
-            } else {
+            }
+            break;
+          case ClpSimplex::atUpperBound:
+            if (value>tolerance) {
               // store square in list
               if (infeas[iSequence+addSequence])
@@ -541 +462 @@
               else
                 infeasible_->quickAdd(iSequence+addSequence,value*value);
+            } else {
+              infeasible_->zero(iSequence+addSequence);
+            }
+            break;
+          case ClpSimplex::atLowerBound:
+            if (value<-tolerance) {
+              // store square in list
+              if (infeas[iSequence+addSequence])
+                infeas[iSequence+addSequence] = value*value; // already there
+              else
+                infeasible_->quickAdd(iSequence+addSequence,value*value);
+            } else {
+              infeasible_->zero(iSequence+addSequence);
             }
           }
         }
-      }
-    }
-  }
-#endif
+      } else {
+        // we can go both ways
+        ClpNonLinearCost * nonLinear = model_->nonLinearCost(); 
+        for (iSequence=0;iSequence<number;iSequence++) {
+          double value = reducedCost[iSequence];
+          ClpSimplex::Status status = model_->getStatus(iSequence+addSequence);
+          
+          switch(status) {
+            
+          case ClpSimplex::basic:
+            infeasible_->zero(iSequence+addSequence);
+          case ClpSimplex::isFixed:
+            break;
+          case ClpSimplex::isFree:
+          case ClpSimplex::superBasic:
+            if (fabs(value)>tolerance) {
+              // we are going to bias towards free (but only if reasonable)
+              value *= FREE_BIAS;
+              // store square in list
+              if (infeas[iSequence+addSequence])
+                infeas[iSequence+addSequence] = value*value; // already there
+              else
+                infeasible_->quickAdd(iSequence+addSequence,value*value);
+            } else {
+              infeasible_->zero(iSequence+addSequence);
+            }
+            break;
+          case ClpSimplex::atUpperBound:
+            if (value>tolerance) {
+              // store square in list
+              if (infeas[iSequence+addSequence])
+                infeas[iSequence+addSequence] = value*value; // already there
+              else
+                infeasible_->quickAdd(iSequence+addSequence,value*value);
+            } else {
+              // look other way - change up should be negative
+              value -= nonLinear->changeUpInCost(iSequence+addSequence);
+              if (value>-tolerance) {
+                infeasible_->zero(iSequence+addSequence);
+              } else {
+                // store square in list
+                if (infeas[iSequence+addSequence])
+                  infeas[iSequence+addSequence] = value*value; // already there
+                else
+                  infeasible_->quickAdd(iSequence+addSequence,value*value);
+              }
+            }
+            break;
+          case ClpSimplex::atLowerBound:
+            if (value<-tolerance) {
+              // store square in list
+              if (infeas[iSequence+addSequence])
+                infeas[iSequence+addSequence] = value*value; // already there
+              else
+                infeasible_->quickAdd(iSequence+addSequence,value*value);
+            } else {
+              // look other way - change down should be positive
+              value -= nonLinear->changeDownInCost(iSequence+addSequence);
+              if (value<tolerance) {
+                infeasible_->zero(iSequence+addSequence);
+              } else {
+                // store square in list
+                if (infeas[iSequence+addSequence])
+                  infeas[iSequence+addSequence] = value*value; // already there
+                else
+                  infeasible_->quickAdd(iSequence+addSequence,value*value);
+              }
+            }
+          }
+        }
+      }
+    }
+  }
   // for weights update we use pivotSequence
   pivotRow = pivotSequence_;
@@ -704 +707 @@
     //weight=1.0;
     if (value>bestDj*weight&&value>tolerance) {
-      // check flagged variable
+      // check flagged variable and correct dj
       if (!model_->flagged(iSequence)) {
         /*if (model_->numberIterations()%100==0)
@@ -1106 +1109 @@
       int nWords = (number+31)>>5;
       reference_ = new unsigned int[nWords];
-      assert (sizeof(unsigned int)==4);
       // tiny overhead to zero out (stops valgrind complaining)
       memset(reference_,0,nWords*sizeof(int));

branches/pre/ClpPrimalQuadraticDantzig.cpp

@@ -90 +90 @@
   //double tolerance=model_->currentPrimalTolerance();
 
-  double bestDj = model_->dualTolerance();
-  double bestInfeasibleDj = 10.0*model_->dualTolerance();
+  double dualTolerance = model_->dualTolerance()*1000.0;
+  double bestDj = dualTolerance;
+  double bestInfeasibleDj = 10.0*dualTolerance;
   int bestSequence=-1;
   double dualIn=0.0;
@@ -105 +106 @@
   const int * columnLength = matrix->getVectorLengths();
   const int * which = quadraticInfo_->quadraticSequence();
-  const double * objective = quadraticInfo_->originalObjective();
+  const double * objective = quadraticInfo_->linearObjective();
+  const double * djWeight = quadraticInfo_->djWeight();
   for (iSequence=0;iSequence<originalNumberColumns;iSequence++) {
     if (model_->flagged(iSequence))
@@ -122 +124 @@
       }
     }
+    double value2 = value*djWeight[iSequence];
+    if (fabs(value2)>dualTolerance)
+      value=value2;
+    else if (value<-dualTolerance)
+      value = -1.001*dualTolerance;
+    else if (value>dualTolerance)
+      value = 1.001*dualTolerance;
+    if (djWeight[iSequence]<1.0e-6)
+      value=value2;
     ClpSimplex::Status status = model_->getStatus(iSequence);
     
@@ -133 +144 @@
         bestSequence = jSequence;
         dualIn = value;
+        abort();
       }
       break;
@@ -171 +183 @@
     // for L rows pi negative okay so choose if positive
     double value = solution[iPi];
+    double value2 = value*djWeight[iSequence];
+    if (fabs(value2)>dualTolerance)
+      value=value2;
+    else if (value<-dualTolerance)
+      value = -1.001*dualTolerance;
+    else if (value>dualTolerance)
+      value = 1.001*dualTolerance;
+    if (djWeight[iSequence]<1.0e-6)
+      value=value2;
     ClpSimplex::Status status = model_->getStatus(iSequence);
     switch(status) {
@@ -208 +229 @@
   if (bestSequence>=0) {
     model_->djRegion()[bestSequence] = dualIn;
-#if 0
-    // free up depending which bound
-    if (solution[iSequence]<trueLower[iSequence]+tolerance) {
-      model_->setColumnStatus(iSequence,ClpSimplex::atLowerBound);
-      upper[iSequence]=trueUpper[iSequence];
-    } else if (solution[iSequence]>trueUpper[iSequence]-tolerance) {
-      model_->setColumnStatus(iSequence,ClpSimplex::atUpperBound);
-      lower[iSequence]=trueLower[iSequence];
-    } else {
-      abort();
-    }
-#endif
   }
   return bestSequence;

branches/pre/ClpSimplex.cpp

@@ -2681 +2681 @@
 int ClpSimplex::strongBranching(int numberVariables,const int * variables,
                                 double * newLower, double * newUpper,
+                                double ** outputSolution,
+                                int * outputStatus, int * outputIterations,
                                 bool stopOnFirstInfeasible,
                                 bool alwaysFinish)
 {
   return ((ClpSimplexDual *) this)->strongBranching(numberVariables,variables,
-                                                    newLower,  newUpper,
+                                                    newLower,  newUpper,outputSolution,
+                                                    outputStatus, outputIterations,
                                                     stopOnFirstInfeasible,
                                                     alwaysFinish);
@@ -4484 +4487 @@
 static bool equalDouble(double value1, double value2)
 {
-  assert(sizeof(unsigned int)*2==sizeof(double));
   unsigned int *i1 = (unsigned int *) &value1;
   unsigned int *i2 = (unsigned int *) &value2;
-  return (i1[0]==i2[0]&&i1[1]==i2[1]);
+  if (sizeof(unsigned int)*2==sizeof(double)) 
+    return (i1[0]==i2[0]&&i1[1]==i2[1]);
+  else
+    return (i1[0]==i2[0]);
 }
 int
@@ -4580 +4585 @@
       return -2;
     } else {
-      model_->messageHandler()->message(CLP_LOOP,model_->messages())
-        <<CoinMessageEol;
+      // look at solution and maybe declare victory
+      if (infeasibility<1.0e-4) {
+        return 0;
+      } else {
+        model_->messageHandler()->message(CLP_LOOP,model_->messages())
+          <<CoinMessageEol;
 #ifndef NDEBUG
-      abort();
+        abort();
 #endif
-      // look at solution and maybe declare victory
-      if (infeasibility<1.0e-4)
-        return 0;
-      else
         return 3;
+      }
     }
   }

branches/pre/ClpSimplexDual.cpp

@@ -2943 +2943 @@
 int ClpSimplexDual::strongBranching(int numberVariables,const int * variables,
                                     double * newLower, double * newUpper,
+                                    double ** outputSolution,
+                                    int * outputStatus, int * outputIterations,
                                     bool stopOnFirstInfeasible,
                                     bool alwaysFinish)
@@ -2969 +2971 @@
   factorization_->slackValue(-1.0);
   factorization_->zeroTolerance(1.0e-13);
-  // save if sparse factorization wanted
-  int saveSparse = factorization_->sparseThreshold();
 
   int factorizationStatus = internalFactorize(0);
@@ -2979 +2979 @@
       <<factorizationStatus
       <<CoinMessageEol;
-  if (saveSparse) {
-    // use default at present
-    factorization_->sparseThreshold(0);
-    factorization_->goSparse();
-  }
   
   // save stuff
@@ -3008 +3003 @@
   // need to save/restore weights.
 
+  int iSolution = 0;
   for (i=0;i<numberVariables;i++) {
     int iColumn = variables[i];
@@ -3024 +3020 @@
     // Start of fast iterations
     int status = fastDual(alwaysFinish);
-
+    if (status) {
+      // not finished - might be optimal
+      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
+      double limit = 0.0;
+      getDblParam(ClpDualObjectiveLimit, limit);
+      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
+          optimizationDirection_*limit) { 
+        problemStatus_=0;
+      } 
+      status=problemStatus_;
+    }
+
+    if (scalingFlag_<=0) {
+      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
+    } else {
+      int j;
+      double * sol = outputSolution[iSolution];
+      for (j=0;j<numberColumns_;j++) 
+        sol[j] = solution_[j]*columnScale_[j];
+    }
+    outputStatus[iSolution]=status;
+    outputIterations[iSolution]=numberIterations_;
+    iSolution++;
     // restore
     memcpy(solution_,saveSolution,
@@ -3061 +3079 @@
     // Start of fast iterations
     status = fastDual(alwaysFinish);
+    if (status) {
+      // not finished - might be optimal
+      checkPrimalSolution(rowActivityWork_,columnActivityWork_);
+      double limit = 0.0;
+      getDblParam(ClpDualObjectiveLimit, limit);
+      if (!numberPrimalInfeasibilities_&&objectiveValue()*optimizationDirection_<
+          optimizationDirection_*limit) { 
+        problemStatus_=0;
+      } 
+      status=problemStatus_;
+    }
+    if (scalingFlag_<=0) {
+      memcpy(outputSolution[iSolution],solution_,numberColumns_*sizeof(double));
+    } else {
+      int j;
+      double * sol = outputSolution[iSolution];
+      for (j=0;j<numberColumns_;j++) 
+        sol[j] = solution_[j]*columnScale_[j];
+    }
+    outputStatus[iSolution]=status;
+    outputIterations[iSolution]=numberIterations_;
+    iSolution++;
 
     // restore
@@ -3121 +3161 @@
   delete [] savePivot;
 
-  factorization_->sparseThreshold(saveSparse);
   // Get rid of some arrays and empty factorization
   deleteRim();
@@ -3138 +3177 @@
   }
   memcpy(saveStatus,status_,(numberColumns_+numberRows_)*sizeof(char));
+  int iSolution =0;
   for (i=0;i<numberVariables;i++) {
     int iColumn = variables[i];
@@ -3146 +3186 @@
     double saveUpper = columnUpper_[iColumn];
     columnUpper_[iColumn] = newUpper[i];
-    dual();
+    int status=dual(0);
+    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
+    outputStatus[iSolution]=status;
+    outputIterations[iSolution]=numberIterations_;
+    iSolution++;
 
     // restore
@@ -3168 +3212 @@
     double saveLower = columnLower_[iColumn];
     columnLower_[iColumn] = newLower[i];
-    dual();
+    status=dual(0);
+    memcpy(outputSolution[iSolution],columnActivity_,numberColumns_*sizeof(double));
+    outputStatus[iSolution]=status;
+    outputIterations[iSolution]=numberIterations_;
+    iSolution++;
 
     // restore
@@ -3224 +3272 @@
 {
   algorithm_ = -1;
+  // save data
+  ClpDataSave data = saveData();
   dualTolerance_=dblParam_[ClpDualTolerance];
   primalTolerance_=dblParam_[ClpPrimalTolerance];
@@ -3239 +3289 @@
   problemStatus_ = -1;
   numberIterations_=0;
+  factorization_->sparseThreshold(0);
+  factorization_->goSparse();
 
   int lastCleaned=0; // last time objective or bounds cleaned up
@@ -3270 +3322 @@
   int returnCode = 0;
 
+  int iRow,iColumn;
   while (problemStatus_<0) {
-    int iRow,iColumn;
     // clear
     for (iRow=0;iRow<4;iRow++) {
@@ -3295 +3347 @@
       double * givenPi=NULL;
       returnCode = whileIterating(givenPi);
-      if (!alwaysFinish&&returnCode<1) {
+      if (!alwaysFinish&&(returnCode<1||returnCode==3)) {
         double limit = 0.0;
         getDblParam(ClpDualObjectiveLimit, limit);
@@ -3301 +3353 @@
            optimizationDirection_*limit|| 
            numberAtFakeBound()) {
+          returnCode=1;
           // can't say anything interesting - might as well return
 #ifdef CLP_DEBUG
@@ -3306 +3359 @@
                  numberIterations_,returnCode);
 #endif
-          returnCode=1;
           break;
         }
@@ -3314 +3366 @@
   }
 
+  // clear
+  for (iRow=0;iRow<4;iRow++) {
+    rowArray_[iRow]->clear();
+  }    
+  
+  for (iColumn=0;iColumn<2;iColumn++) {
+    columnArray_[iColumn]->clear();
+  }    
   assert(!numberFake_||returnCode||problemStatus_); // all bounds should be okay
+  // Restore any saved stuff
+  restoreData(data);
   dualBound_ = saveDualBound;
   return returnCode;

branches/pre/ClpSimplexPrimalQuadratic.cpp

@@ -9 +9 @@
 #include "ClpSimplexPrimalQuadratic.hpp"
 #include "ClpPrimalQuadraticDantzig.hpp"
+#include "ClpPrimalColumnDantzig.hpp"
+#include "ClpQuadraticObjective.hpp"
 #include "ClpFactorization.hpp"
 #include "ClpNonLinearCost.hpp"
@@ -73 +75 @@
 
   // Get list of non linear columns
-  CoinPackedMatrix * quadratic = quadraticObjective();
+  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
+  CoinPackedMatrix * quadratic = NULL;
+  if (quadraticObj)
+    quadratic = quadraticObj->quadraticObjective();
   if (!quadratic) {
     // no quadratic part
@@ -496 +501 @@
 #endif  
   // Now do quadratic
-  ClpPrimalQuadraticDantzig dantzigP(model2,&info);
-  model2->setPrimalColumnPivotAlgorithm(dantzigP);
+  //ClpPrimalQuadraticDantzig dantzigP(model2,&info);
+  //ClpPrimalColumnDantzig dantzigP;
+  //model2->setPrimalColumnPivotAlgorithm(dantzigP);
   model2->messageHandler()->setLogLevel(63);
   model2->primalQuadratic2(&info,phase);
@@ -520 +526 @@
   const int * which = info->quadraticSequence();
   //const int * back = info->backSequence();
-  assert (!scalingFlag_);
   // initialize - values pass coding and algorithm_ is +1
   if (!startup(1)) {
@@ -529 +534 @@
     info->setSequenceIn(-1);
     info->setCrucialSj(-1);
+    bool deleteCosts=false;
+    if (scalingFlag_) {
+      // get own copy of original objective and scale
+      ClpQuadraticObjective * quadraticObj = new ClpQuadraticObjective(*info->originalObjective());
+      info->setOriginalObjective(quadraticObj);
+      CoinPackedMatrix * quadratic = info->quadraticObjective();
+      double * objective = info->linearObjective();
+      // replace column by column
+      double * newElement = new double[numberXColumns];
+      // scale column copy
+      // get matrix data pointers
+      const int * row = quadratic->getIndices();
+      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
+      const int * columnLength = quadratic->getVectorLengths(); 
+      const double * elementByColumn = quadratic->getElements();
+      double direction = optimizationDirection_;
+      // direction is actually scale out not scale in
+      if (direction)
+        direction = 1.0/direction;
+      int iColumn;
+      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
+        int j;
+        double scale = columnScale_[iColumn]*direction;
+        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
+        const int * columnsInThisColumn = row + columnStart[iColumn];
+        int number = columnLength[iColumn];
+        assert (number<=numberXColumns);
+        for (j=0;j<number;j++) {
+          int iColumn2 = columnsInThisColumn[j];
+          newElement[j] = elementsInThisColumn[j]*scale*columnScale_[iColumn2];
+        }
+        quadratic->replaceVector(iColumn,number,newElement);
+        objective[iColumn] *= direction*columnScale_[iColumn];
+      }
+      delete [] newElement;
+      deleteCosts=true;
+    } else if (optimizationDirection_!=1.0) {
+      // get own copy of original objective and scale
+      ClpQuadraticObjective * quadraticObj = new ClpQuadraticObjective(*info->originalObjective());
+      info->setOriginalObjective(quadraticObj);
+      CoinPackedMatrix * quadratic = info->quadraticObjective();
+      double * objective = info->linearObjective();
+      // replace column by column
+      double * newElement = new double[numberXColumns];
+      // get matrix data pointers
+      const CoinBigIndex * columnStart = quadratic->getVectorStarts();
+      const int * columnLength = quadratic->getVectorLengths(); 
+      const double * elementByColumn = quadratic->getElements();
+      double direction = optimizationDirection_;
+      // direction is actually scale out not scale in
+      if (direction)
+        direction = 1.0/direction;
+      int iColumn;
+      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
+        int j;
+        const double * elementsInThisColumn = elementByColumn + columnStart[iColumn];
+        int number = columnLength[iColumn];
+        assert (number<=numberXColumns);
+        for (j=0;j<number;j++) {
+          newElement[j] = elementsInThisColumn[j]*direction;
+        }
+        quadratic->replaceVector(iColumn,number,newElement);
+        objective[iColumn] *= direction;
+      }
+      delete [] newElement;
+      deleteCosts=true;
+    }
     
     // Say no pivot has occurred (for steepest edge and updates)
@@ -568 +640 @@
         break; // declare victory
       
-      // 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();
-      objectiveValue_=0.0;
-      for (iColumn=0;iColumn<numberXColumns;iColumn++) {
-        double valueI = solution_[iColumn];
-        objectiveValue_ += valueI*originalCost[iColumn];
-        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;
-        }
-      }
+      checkComplementarity (info,rowArray_[0],rowArray_[1]);
 
       // Say good factorization
@@ -656 +708 @@
 
       // exit if victory declared
+      dualIn_=0.0; // so no updates
       if (!info->currentPhase()&&info->sequenceIn()<0&&primalColumnPivot_->pivotColumn(rowArray_[1],
                                           rowArray_[2],rowArray_[3],
@@ -667 +720 @@
       problemStatus_=-1;
       whileIterating(info);
+    }
+    if (deleteCosts) {
+      // delete scaled copy of objective
+      delete info->originalObjective();
     }
   }
@@ -697 +754 @@
   double saveObjective = objectiveValue_;
   int numberXColumns = info->numberXColumns();
-  int numberXRows = info->numberXRows();
   const int * which = info->quadraticSequence();
   const double * pi = solution_+numberXColumns;
-  // Matrix for linear stuff
-  const int * row = matrix_->getIndices();
-  const int * columnStart = matrix_->getVectorStarts();
-  const double * element =  matrix_->getElements();
-  const int * columnLength = matrix_->getVectorLengths();
   int nextSequenceIn=info->sequenceIn();
   int oldSequenceIn=nextSequenceIn;
@@ -771 +822 @@
       } else {
         saveSequenceIn=sequenceIn_;
+        createDjs(info,rowArray_[1],rowArray_[2]);
+        dualIn_=0.0; // so no updates
         primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
                      columnArray_[0],columnArray_[1]);
@@ -777 +830 @@
       saveSequenceIn=sequenceIn_;
       sequenceIn_ = nextSequenceIn;
-      if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_) 
-        cleanupIteration=false;
-      else
+      if (phase==2) {
+        if (sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_) 
+          cleanupIteration=false;
+        else
+          cleanupIteration=true;
+      } else {
         cleanupIteration=true;
+      }
     }
     pivotRow_=-1;
@@ -786 +843 @@
     rowArray_[1]->clear();
     if (sequenceIn_>=0) {
-      if (numberIterations_>=8796000)
-        printf("loxx 4721 %g %g\n",lower_[4721],upper_[4721]);
       nextSequenceIn=-1;
       // we found a pivot column
@@ -793 +848 @@
       // 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);
-        }
+        checkComplementarity (info,rowArray_[3],rowArray_[1]);
+        printf("True objective is %g, infeas cost %g, sum %g\n",
+               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
         objectiveValue_=saveObjective;
         returnCode=-1;
@@ -859 +861 @@
         chosen=-1;
         unpack(rowArray_[1]);
+        // compute dj in case linear
+        {
+          dualIn_=cost_[sequenceIn_];
+          int j;
+          const double * element = rowArray_[1]->denseVector();
+          int numberNonZero=rowArray_[1]->getNumElements();
+          const int * whichRow = rowArray_[1]->getIndices();
+          bool packed = rowArray_[1]->packedMode();
+          if (packed) {
+            for (j=0;j<numberNonZero; j++) {
+              int iRow = whichRow[j];
+              dualIn_ -= element[j]*pi[iRow];
+            }
+          } else {
+            for (j=0;j<numberNonZero; j++) {
+              int iRow = whichRow[j];
+              dualIn_ -= element[iRow]*pi[iRow];
+            }
+          }
+        }
         factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
         if (cleanupIteration) {
@@ -918 +940 @@
             } 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];
-              }
+              //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 {
@@ -957 +979 @@
                              info,
                              cleanupIteration);
-        if (pivotRow_==-1&&phase==2&&fabs(dualIn_)<1.0e-3) {
+        if (pivotRow_==-1&&(phase==2||cleanupIteration)&&fabs(dualIn_)<1.0e-3) {
           // try other way
           dualIn_=-dualIn_;
@@ -963 +985 @@
           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],
@@ -1139 +1158 @@
           printf("%d Sj value went from %g to %g\n",crucialSj,oldSj,solution_[info->crucialSj()]);
         }
-        {
-          double oldValue = objectiveValue_;
-          // 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;
-            }
+        double oldObjectiveValue = objectiveValue_;
+        checkComplementarity (info,rowArray_[2],rowArray_[3]);
+        printf("After Housekeeping True objective is %g, infeas cost %g, sum %g\n",
+               objectiveValue_,info->infeasCost(),objectiveValue_+info->infeasCost());
+        if (objectiveValue_>oldObjectiveValue) {
+          // make less likely this one will come in again
+          double multiplier = CoinDrand48();
+          int iSequence;
+          if (!cleanupIteration) {
+            iSequence = sequenceIn_;
+          } else {
+            iSequence = info->crucialSj();
+            if (iSequence>numberRows_)
+              iSequence -= numberRows_;
+            else
+              iSequence = numberColumns_ +(iSequence-numberXColumns);
           }
-          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 (objectiveValue_>oldValue) {
-            // make less likely this one will come in again
-            double multiplier = CoinDrand48();
-            int iSequence;
-            if (!cleanupIteration) {
-              iSequence = sequenceIn_;
-            } else {
-              int iSequence = info->crucialSj();
-              if (iSequence>numberRows_)
-                iSequence -= numberRows_;
-              else
-                iSequence = numberColumns_ +(iSequence-numberXColumns);
-            }
-            info->djWeight()[iSequence] /= (5.0+multiplier);
-          }
+          info->djWeight()[iSequence] /= (5.0+multiplier);
         }
         if (sequenceOut_>=numberXColumns&&sequenceOut_<numberColumns_) {
@@ -1217 +1183 @@
             info->setCrucialSj(-1);
         }
-        checkComplementarity (info,rowArray_[2],rowArray_[3]);
 
         if (whatNext==1) {
@@ -1233 +1198 @@
         // may not be correct on second time
         if (result&&(returnCode==-1||returnCode==-2||returnCode==-3)) {
-          assert(sequenceOut_<numberXColumns||
-                 sequenceOut_>=numberColumns_);
           int crucialSj=info->crucialSj();
           if (sequenceOut_<numberXColumns) {
             chosen =sequenceOut_ + numberRows_; // sj variable
-          } else {
+          } else if (sequenceOut_>=numberColumns_) {
             // Does this mean we can change pi
             int iRow = sequenceOut_-numberColumns_;
@@ -1250 +1213 @@
               abort();
             }
+          } else if (sequenceOut_<numberRows_) {
+            // pi out
+            chosen = sequenceOut_-numberXColumns+numberColumns_;
+          } else {
+            // sj out
+            chosen = sequenceOut_-numberRows_;
           }
           // But double check as original incoming might have gone out
@@ -1381 +1350 @@
   // Work out coefficients for quadratic term
   // This is expanded one
-  const CoinPackedMatrix * quadratic = quadraticObjective();
+  const CoinPackedMatrix * quadratic = info->quadraticObjective();
   const int * columnQuadratic = quadratic->getIndices();
   const int * columnQuadraticStart = quadratic->getVectorStarts();
@@ -1420 +1389 @@
       //printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
     } else {
-      coeffSlack += alpha*cost_[iPivot];
-      //printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
-      if (iPivot==crucialSj) {
-        tj = alpha;
-        iSjRow = iRow;
+      if (iPivot>=numberColumns_) {
+        // ? do we go through column of pi
+        assert (iPivot!=crucialSj);
+        coeffSlack += alpha*cost_[iPivot];
+      } else if (iPivot<numberRows_) {
+        // ? do we go through column of pi
+        if (iPivot==crucialSj) {
+          tj = alpha;
+          iSjRow = iRow;
+        }
+      } else {
+        if (iPivot==crucialSj) {
+          tj = alpha;
+          iSjRow = iRow;
+        }
       }
     }
@@ -1441 +1420 @@
   }
   rhsArray->setNumElements(number2);
-  if (numberIterations_!=-701) {  //if (numberIterations_>=0||cleanupIteration) {
-    for (iIndex=0;iIndex<number2;iIndex++) {
-      int iColumn=index2[iIndex];
-      //double valueI = solution_[iColumn];
-      double alphaI = rhs[iColumn];
-      int j;
-      for (j=columnQuadraticStart[iColumn];
-           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
-        int jColumn = columnQuadratic[j];
-        double valueJ = solution_[jColumn];
-        double alphaJ = rhs[jColumn];
-        double elementValue = quadraticElement[j];
-        coeff1 += (valueJ*alphaI)*elementValue;
-        coeff2 += (alphaI*alphaJ)*elementValue;
-      }
-    }
-  } else {
-    const int * row = matrix_->getIndices();
-    const int * columnStart = matrix_->getVectorStarts();
-    const int * columnLength = matrix_->getVectorLengths();
-    const double * element = matrix_->getElements();
+  for (iIndex=0;iIndex<number2;iIndex++) {
+    int iColumn=index2[iIndex];
+    //double valueI = solution_[iColumn];
+    double alphaI = rhs[iColumn];
     int j;
-    int jRow = sequenceIn_+info->numberXRows();
-    printf("sequence in %d, cost %g\n",sequenceIn_,
-           upper_[jRow+numberColumns_]);
-    double xx=0.0;
-    for (j=0;j<numberColumns_;j++) {
-      int k;
-      for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
-        if (row[k]==jRow) {
-          printf ("col %d, el %g, sol %g, contr %g\n",
-                  j,element[k],solution_[j],element[k]*solution_[j]);
-          xx+= element[k]*solution_[j];
-        }
-      }
-    }
-    printf("sum %g\n",xx);
-    for (iIndex=0;iIndex<number2;iIndex++) {
-      int iColumn=index2[iIndex];
-      double valueI = solution_[iColumn];
-      double alphaI = rhs[iColumn];
-      //rhs[iColumn]=0.0;
-      printf("Column %d, alpha %g sol %g cost %g, contr %g\n",
-             iColumn,alphaI,valueI,cost_[iColumn],
-             alphaI*cost_[iColumn]);
-      int j;
-      for (j=columnQuadraticStart[iColumn];
-           j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
-        int jColumn = columnQuadratic[j];
-        double valueJ = solution_[jColumn];
-        double alphaJ = rhs[jColumn];
-        double elementValue = quadraticElement[j];
-        printf("j %d alphaJ %g solJ %g el %g, contr %g\n",
-               jColumn,alphaJ,valueJ,elementValue,
-               (valueJ*alphaI)*elementValue);
-        coeff1 += (valueJ*alphaI)*elementValue;
-        coeff2 += (alphaI*alphaJ)*elementValue;
-      }
+    for (j=columnQuadraticStart[iColumn];
+         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
+      int jColumn = columnQuadratic[j];
+      double valueJ = solution_[jColumn];
+      double alphaJ = rhs[jColumn];
+      double elementValue = quadraticElement[j];
+      coeff1 += (valueJ*alphaI)*elementValue;
+      coeff2 += (alphaI*alphaJ)*elementValue;
     }
   }
@@ -1508 +1442 @@
   int accuracyFlag=0;
   if (!cleanupIteration) {
-    //assert (fabs(way*coeff1-dualIn_)<1.0e-4*(1.0+fabs(dualIn_)));
+    if (fabs(dualIn_)<dualTolerance_&&fabs(coeff1)<dualTolerance_) {
+      dualIn_=0.0;
+      coeff1=0.0;
+    }
+    //assert (fabs(way*coeff1-dualIn_)<1.0e-1*(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-4*(1.0+fabs(dualIn_))) {
+    } else if (fabs(way*coeff1-dualIn_)>1.0e-3*(1.0+fabs(dualIn_))) {
       // not wondeful
       accuracyFlag=1;
@@ -1554 +1492 @@
     <<coeff1<<coeff2<<coeffSlack<<dualIn_<<d1<<d2
     <<CoinMessageEol;
-  if (!cleanupIteration)
-    dualIn_ = way*coeff1;
+  //if (!cleanupIteration)
+  //dualIn_ = way*coeff1;
   double mind1d2=1.0e50;
   if (cleanupIteration) {
@@ -1570 +1508 @@
   }
   maximumMovement = min(maximumMovement,mind1d2);
+  double trueMaximumMovement;
+  if (way>0.0) 
+    trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
+  else
+    trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
   rhsArray->clear();
   double tentativeTheta = maximumMovement;
   double upperTheta = maximumMovement;
   bool throwAway=false;
-  if (numberIterations_>=2007) {
+  if (numberIterations_==126) {
     printf("Bad iteration coming up after iteration %d\n",numberIterations_);
   }
 
+  double minimumAny=1.0e50;
+  int whichMinimum=-1;
+  double minimumAlpha=0.0;
   for (iIndex=0;iIndex<number;iIndex++) {
 
@@ -1595 +1541 @@
       oldValue = bound-oldValue;
     }
+    if (iPivot>=numberXColumns&&iPivot<numberColumns_) {
+      double bound=0.0;
+      double oldValue2 = solution(iPivot);
+      if (alpha>0.0) {
+        // basic variable going towards lower bound
+        oldValue2 -= bound;
+      } else if (alpha<0.0) {
+        // basic variable going towards upper bound
+        oldValue2 = bound-oldValue;
+      }
+      double value = oldValue2-minimumAny*fabs(alpha);
+      if (value*oldValue2<0.0) {
+        double ratio = fabs(oldValue2/alpha);
+        if (ratio<minimumAny&&fabs(alpha)>1.0e2*acceptablePivot) {
+          minimumAny = ratio;
+          whichMinimum = iRow;
+          minimumAlpha= fabs(alpha);
+        }
+      }
+    }
     double value = oldValue-tentativeTheta*fabs(alpha);
     assert (oldValue>=-primalTolerance_*1.0001);
@@ -1616 +1582 @@
 
   bool goBackOne = false;
+  double fake=1.0e-2;
 
   if (numberRemaining) {
@@ -1669 +1636 @@
       //dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
       //dualCheck=0.0;
-      if (totalThru+thruThis>=dualCheck) {
+      if ((totalThru+thruThis>=dualCheck&&bestPivot>acceptablePivot)||fake*bestPivot>1.0e-3) {
         // We should be pivoting in this batch
         // so compress down to this lot
@@ -1748 +1715 @@
           } else {
             dualCheck = - 2.0*coeff2*theta_ - coeff1;
-            //dualCheck =0.0;
-            if (totalThru>=dualCheck&&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
+            if ((totalThru>=dualCheck||fake*bestPivot>1.0e-3)
+                &&(sequenceIn_<numberXColumns||sequenceIn_>=numberColumns_)) {
               if (!cleanupIteration) {
                 // We can pivot out sj
                 if (upperTheta==maximumMovement) {
                   printf("figures\n");
-                } else {
-                  if (lastThru>=dualCheck)
+                } else if (iSjRow2>=0) {
+                  if (lastThru>=dualCheck&&theta_<trueMaximumMovement) {
                     printf("totalThru %g, lastThru %g - taking sj to zero?\n",totalThru,lastThru);
+                    // make sure will take correct path
+                    mind1d2=1.0;
+                    maximumMovement=1.0;
+                    d2=0.0;
+                    pivotRow_=-1;
+                  }
                 }
               } else {
@@ -1763 +1736 @@
               }
               break; // no point trying another loop
-            } else if (totalThru>=dualCheck||upperTheta==maximumMovement) {
+            } else if (totalThru>=dualCheck||fake*bestPivot>1.0e-3||upperTheta==maximumMovement) {
               break; // no point trying another loop
             } else {
@@ -1807 +1780 @@
     }
   }
+  if (0&&minimumAny<theta_&&cleanupIteration&&bestEverPivot<1.0e-2
+      &&minimumAlpha>bestEverPivot*1.0e2&&minimumAny>1.0e-1&&info->currentPhase()==0) {
+    printf("Possible other pivot %d %g %g - alpha %g\n",
+           whichMinimum,minimumAny,theta_,minimumAlpha);
+    pivotRow_ = whichMinimum;
+    alpha_ = work[pivotRow_];
+    // translate to sequence
+    sequenceOut_ = pivotVariable_[pivotRow_];
+    valueOut_ = solution(sequenceOut_);
+    lowerOut_=0.0;
+    upperOut_=0.0;
+    theta_= fabs(valueOut_/alpha_);
+
+    if (way<0.0) 
+      theta_ = - theta_;
+    if (alpha_*way<0.0) {
+      directionOut_=-1;      // to upper bound
+    } else {
+      directionOut_=1;      // to lower bound
+    }
+    dualOut_ = reducedCost(sequenceOut_);
+  } else {
 
   if (pivotRow_>=0) {
-    
+    if (0) {
+      double move = coeff1*theta_+coeff2*theta_*theta_;
+      printf("Predicted change in obj is %g %g %g\n",
+             move,objectiveValue_-move,objectiveValue_+move);
+    }
 #define MINIMUMTHETA 1.0e-12
     // will we need to increase tolerance
@@ -1865 +1864 @@
       accuracyFlag=2;
     }
-    double trueMaximumMovement;
-    if (way>0.0) 
-      trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
-    else
-      trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
     if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
       // flip
@@ -1898 +1892 @@
       } else {
         // incoming dj going to zero
-        assert(!cleanupIteration);
-        if (!cleanupIteration)
+        if (!cleanupIteration) {
           result=0;
-        iSjRow=iSjRow2;
-        crucialSj = pivotVariable_[iSjRow];
-      }
-      assert (pivotRow_<0);
-      setColumnStatus(crucialSj,isFree);
-      pivotRow_ = iSjRow;
-      alpha_ = work[pivotRow_];
-      // translate to sequence
-      sequenceOut_ = pivotVariable_[pivotRow_];
-      assert (sequenceOut_==crucialSj);
-      valueOut_ = solution(sequenceOut_);
-      theta_ = maximumMovement;
-      if (way<0.0) 
-        theta_ = - theta_;
-      lowerOut_=0.0;
-      upperOut_=0.0;
-      //????
-      dualOut_ = reducedCost(sequenceOut_);
+          iSjRow=iSjRow2;
+          crucialSj = pivotVariable_[iSjRow];
+          assert (pivotRow_<0);
+        } else {
+          assert (fabs(dualIn_)<1.0e-3);
+          result=1;
+          if (minimumAlpha>0.0) {
+            // could do this in other places if pivot looks good
+            printf("Should take minimum\n");
+            pivotRow_ = whichMinimum;
+            alpha_ = work[pivotRow_];
+            // translate to sequence
+            sequenceOut_ = pivotVariable_[pivotRow_];
+            valueOut_ = solution(sequenceOut_);
+            theta_ = minimumAny;
+            if (way<0.0) 
+              theta_ = - theta_;
+            lowerOut_=0.0;
+            upperOut_=0.0;
+            //????
+            dualOut_ = reducedCost(sequenceOut_);
+          }
+        }
+      }
+      if (!result) {
+        setColumnStatus(crucialSj,isFree);
+        pivotRow_ = iSjRow;
+        alpha_ = work[pivotRow_];
+        // translate to sequence
+        sequenceOut_ = pivotVariable_[pivotRow_];
+        assert (sequenceOut_==crucialSj);
+        valueOut_ = solution(sequenceOut_);
+        theta_ = fabs(valueOut_/alpha_);
+        assert (fabs(maximumMovement-theta_)<1.0e-3*1.0+maximumMovement);
+        if (way<0.0) 
+          theta_ = - theta_;
+        lowerOut_=0.0;
+        upperOut_=0.0;
+        //????
+        dualOut_ = reducedCost(sequenceOut_);
+      }
     } else {
       // need to do something
@@ -1924 +1940 @@
     }
   }
+  }
+
 
   // clear arrays
@@ -1940 +1958 @@
     objectiveValue_ =0.0;
     CoinPackedMatrix * quadratic = 
-      info->originalModel()->quadraticObjective();
+      info->originalObjective()->quadraticObjective();
     const int * columnQuadratic = quadratic->getIndices();
     const int * columnQuadraticStart = quadratic->getVectorStarts();
@@ -1946 +1964 @@
     const double * quadraticElement = quadratic->getElements();
     int numberColumns = info->numberXColumns();
-    const double * objective = info->originalObjective();
+    const double * objective = info->linearObjective();
     for (iColumn=0;iColumn<numberColumns;iColumn++) 
       objectiveValue_ += objective[iColumn]*solution_[iColumn];
@@ -2073 +2091 @@
   }
   // Check if looping
-  // Take out for now
   int loop;
-  if (type!=2&&0) 
-    loop = progress->looping();
+  if (type!=2) 
+    loop = progress->looping(); // saves iteration number
   else
     loop=-1;
+  //loop=-1;
   if (loop>=0) {
     problemStatus_ = loop; //exit if in loop
@@ -2107 +2125 @@
   // give code benefit of doubt
   if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
-      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
+      sumOfRelaxedPrimalInfeasibilities_ == 0.0&&info->sequenceIn()<0) {
     // say optimal (with these bounds etc)
     numberDualInfeasibilities_ = 0;
@@ -2332 +2350 @@
 {
   const int * which = info->quadraticSequence();
-  // Matrix for linear stuff
-  const int * row = matrix_->getIndices();
-  const int * columnStart = matrix_->getVectorStarts();
-  const double * element =  matrix_->getElements();
-  const int * columnLength = matrix_->getVectorLengths();
   int numberXColumns = info->numberXColumns();
   int numberXRows = info->numberXRows();
@@ -2343 +2356 @@
   int start=numberXColumns+numberXRows;
   int jSequence;
+  const double * djWeight = info->djWeight();
   // Actually only need to do this after invert (use extra parameter to createDjs)
   // or when infeasibilities change
@@ -2458 +2472 @@
     }
   }
+  // fill in linear ones
+  memcpy(dj_,cost_,numberXColumns*sizeof(double));
+  if (!rowScale_) {
+    matrix_->transposeTimes(-1.0,pi,dj_);
+  } else {
+    matrix_->transposeTimes(-1.0,pi,dj_,rowScale_,columnScale_);
+  }
+  memset(dj_+numberXColumns,0,(numberXRows+info->numberQuadraticColumns())*sizeof(double));
   for (iSequence=0;iSequence<numberXColumns;iSequence++) {
     int jSequence = which[iSequence];
@@ -2465 +2487 @@
       value = solution_[jSequence];
     } else {
-      value=0.0;
-      //value=cost_[iSequence];
-      int j;
-      for (j=columnStart[iSequence];j<columnStart[iSequence]+columnLength[iSequence]; j++) {
-        int iRow = row[j];
-        value -= element[j]*pi[iRow];
-      }
-    }
+      value=dj_[iSequence];
+    }
+    double value2 = value*djWeight[iSequence];
+    if (fabs(value2)>dualTolerance_)
+      value=value2;
+    else if (value<-dualTolerance_)
+      value = -1.001*dualTolerance_;
+    else if (value>dualTolerance_)
+      value = 1.001*dualTolerance_;
+    if (djWeight[iSequence]<1.0e-6)
+      value=value2;
     dj_[iSequence]=value;
   }
@@ -2480 +2505 @@
     int iSequence  = jSequence + numberColumns_;
     int iPi = jSequence+numberXColumns;
-    double value = solution_[iPi];
-    dj_[iSequence]=value+cost_[iSequence];
+    double value = solution_[iPi]+cost_[iSequence];
+    double value2 = value*djWeight[iSequence];
+    if (fabs(value2)>dualTolerance_)
+      value=value2;
+    else if (value<-dualTolerance_)
+      value = -1.001*dualTolerance_;
+    else if (value>dualTolerance_)
+      value = 1.001*dualTolerance_;
+    if (djWeight[iSequence]<1.0e-6)
+      value=value2;
+    dj_[iSequence]=value;
   }
 }
 
 int 
-ClpSimplexPrimalQuadratic::checkComplementarity (const  ClpQuadraticInfo * info,
+ClpSimplexPrimalQuadratic::checkComplementarity (ClpQuadraticInfo * info,
                             CoinIndexedVector * array1, CoinIndexedVector * array2)
 {
@@ -2498 +2532 @@
   // Compute objective function from scratch
   const CoinPackedMatrix * quadratic = 
-    info->originalModel()->quadraticObjective();
+    info->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();
+  const double * originalCost = info->linearObjective();
   int iColumn;
   objectiveValue_=0.0;
   double infeasCost=0.0;
   double linearCost=0.0;
+  int number0=0,number1=0,number2=0;
   for (iColumn=0;iColumn<numberXColumns;iColumn++) {
     double valueI = solution_[iColumn];
@@ -2531 +2566 @@
     }
     int jSequence = which[iColumn];
+    if (jSequence>=0) 
+      jSequence += start;
     double value=dj_[iColumn];
     ClpSimplex::Status status = getColumnStatus(iColumn);
+    if (status!=basic&&jSequence>=0) {
+      if (getColumnStatus(jSequence)==basic) 
+        number1++;
+      else
+        number0++;
+    }
     
     switch(status) {
@@ -2538 +2581 @@
     case ClpSimplex::basic:
       if (jSequence>=0) {
-        jSequence += start;
-        if (getColumnStatus(jSequence)==basic)
+        if (getColumnStatus(jSequence)==basic) {
           handler_->message(CLP_QUADRATIC_BOTH,messages_)
             <<"Structural"<<iColumn
             <<solution_[iColumn]<<jSequence<<solution_[jSequence]
             <<CoinMessageEol;
+          number2++;
+        } else {
+          number1++;
+        }
       }
       break;
@@ -2586 +2632 @@
     double value=dj_[iRow+numberColumns_];
     ClpSimplex::Status status = getRowStatus(iRow);
+    if (status!=basic) {
+      if (getColumnStatus(jSequence)==basic) 
+        number1++;
+      else
+        number0++;
+    }
     
     switch(status) {
       
     case ClpSimplex::basic:
-      if (getColumnStatus(jSequence)==basic)
+      if (getColumnStatus(jSequence)==basic) {
         printf("Row %d (%g) and %d (%g) both basic\n",
                iRow,solution_[iColumn],jSequence,solution_[jSequence]);
+          number2++;
+        } else {
+          number1++;
+        }
       break;
     case ClpSimplex::isFixed:
@@ -2616 +2672 @@
     }
   }
+  //printf("number 0 - %d, 1 - %d, 2 - %d\n",number0,number1,number2);
   objectiveValue_ += linearCost+infeasCost;
   assert (infeasCost>=0.0);
   sumOfRelaxedDualInfeasibilities_ =sumDualInfeasibilities_;
+  // But not zero if cleanup iteration
+  if (info->sequenceIn()>=0&&!numberDualInfeasibilities_)
+    numberDualInfeasibilities_=1;
   numberDualInfeasibilitiesWithoutFree_= numberDualInfeasibilities_;
+  info->setInfeasCost(infeasCost);
   return numberDualInfeasibilities_;
 }
@@ -2631 +2692 @@
 
   // Get list of non linear columns
-  CoinPackedMatrix * quadratic = quadraticObjective();
+  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
+  assert (quadraticObj);
+  info.setOriginalObjective(quadraticObj);
+  CoinPackedMatrix * quadratic = info.quadraticObjective();
   if (!quadratic||!quadratic->getNumElements()) {
     // no quadratic part
@@ -2640 +2704 @@
   double * columnLower = this->columnLower();
   double * columnUpper = this->columnUpper();
-  double * objective = this->objective();
+  double * objective = info.linearObjective();
   double * solution = this->primalColumnSolution();
   double * dj = this->dualColumnSolution();
@@ -2982 +3046 @@
   for (iColumn=0;iColumn<numberColumns_;iColumn++) 
     objValue += objective[iColumn]*solution2[iColumn];
-  CoinPackedMatrix * quadratic = quadraticObjective();
+  ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
+  assert (quadraticObj);
+  CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
   double * dj = dualColumnSolution();
   // Matrix for linear stuff
@@ -3053 +3119 @@
 /// Default constructor. 
 ClpQuadraticInfo::ClpQuadraticInfo()
-  : originalModel_(NULL),
+  : originalObjective_(NULL),
     quadraticSequence_(NULL),
     backSequence_(NULL),
@@ -3067 +3133 @@
     numberXRows_(-1),
     numberXColumns_(-1),
-    numberQuadraticColumns_(0)
+    numberQuadraticColumns_(0),
+    infeasCost_(0.0)
 {
 }
 // Constructor from original model
 ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
-  : originalModel_(model),
+  : originalObjective_(NULL),
     quadraticSequence_(NULL),
     backSequence_(NULL),
@@ -3086 +3153 @@
     numberXRows_(-1),
     numberXColumns_(-1),
-    numberQuadraticColumns_(0)
+    numberQuadraticColumns_(0),
+    infeasCost_(0.0)
 {
-  if (originalModel_) {
-    numberXRows_ = originalModel_->numberRows();
-    numberXColumns_ = originalModel_->numberColumns();
+  if (model) {
+    numberXRows_ = model->numberRows();
+    numberXColumns_ = model->numberColumns();
+    originalObjective_ = (dynamic_cast< ClpQuadraticObjective*>(model->objectiveAsObject()));
+    assert (originalObjective_);
     quadraticSequence_ = new int[numberXColumns_];
     backSequence_ = new int[numberXColumns_];
@@ -3118 +3188 @@
 // Copy
 ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
-  : originalModel_(rhs.originalModel_),
+  : originalObjective_(rhs.originalObjective_),
     quadraticSequence_(NULL),
     backSequence_(NULL),
@@ -3129 +3199 @@
     numberXRows_(rhs.numberXRows_),
     numberXColumns_(rhs.numberXColumns_),
-    numberQuadraticColumns_(rhs.numberQuadraticColumns_)
+    numberQuadraticColumns_(rhs.numberQuadraticColumns_),
+    infeasCost_(rhs.infeasCost_)
 {
   if (numberXColumns_) {
@@ -3166 +3237 @@
 {
   if (this != &rhs) {
-    originalModel_ = rhs.originalModel_;
+    originalObjective_ = rhs.originalObjective_;
     delete [] quadraticSequence_;
     delete [] backSequence_;
@@ -3180 +3251 @@
     numberXRows_ = rhs.numberXRows_;
     numberXColumns_ = rhs.numberXColumns_;
+    infeasCost_=rhs.infeasCost_;
     numberQuadraticColumns_=rhs.numberQuadraticColumns_;
     if (numberXColumns_) {
@@ -3258 +3330 @@
   }
 }
+// Quadratic objective
+CoinPackedMatrix * 
+ClpQuadraticInfo::quadraticObjective() const     
+{ 
+  return originalObjective_->quadraticObjective();
+}
+// Linear objective
+double * 
+ClpQuadraticInfo::linearObjective() const     
+{ 
+  return originalObjective_->linearObjective();
+}

branches/pre/Makefile.Clp

@@ -6 +6 @@
 # between then specify the exact level you want, e.g., -O1 or -O2
 OptLevel := -g
-#OptLevel := -O3
+OptLevel := -O3
 
 
@@ -28 +28 @@
 LIBSRC += ClpPrimalColumnPivot.cpp
 LIBSRC += ClpPrimalColumnSteepest.cpp
+LIBSRC += ClpQuadraticObjective.cpp
 LIBSRC += ClpSimplex.cpp
 LIBSRC += ClpSimplexDual.cpp
 LIBSRC += ClpSimplexPrimal.cpp
 LIBSRC += ClpSimplexPrimalQuadratic.cpp
-LIBSRC += ClpPrimalQuadraticDantzig.cpp
+#LIBSRC += ClpPrimalQuadraticDantzig.cpp
 # and Presolve stuff
 LIBSRC += ClpPresolve.cpp

branches/pre/Test/ClpMain.cpp

@@ -1703 +1703 @@
                 double * primalRowSolution = 
                   models[iModel].primalRowSolution();
+                double * rowLower = models[iModel].rowLower();
+                double * rowUpper = models[iModel].rowUpper();
+                double primalTolerance = models[iModel].primalTolerance();
                 char format[6];
                 sprintf(format,"%%-%ds",max(lengthName,8));
                 for (iRow=0;iRow<numberRows;iRow++) {
+                  if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance||
+                      primalRowSolution[iRow]<rowLower[iRow]-primalTolerance)
+                    fprintf(fp,"** ");
                   fprintf(fp,"%7d ",iRow);
                   if (lengthName)
@@ -1718 +1724 @@
                 double * primalColumnSolution = 
                   models[iModel].primalColumnSolution();
+                double * columnLower = models[iModel].columnLower();
+                double * columnUpper = models[iModel].columnUpper();
                 for (iColumn=0;iColumn<numberColumns;iColumn++) {
+                  if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance||
+                      primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance)
+                    fprintf(fp,"** ");
                   fprintf(fp,"%7d ",iColumn);
                   if (lengthName)

branches/pre/Test/unitTest.cpp

@@ -936 +936 @@
   }
   // test network 
-#define QUADRATIC
+  //#define QUADRATIC
 #ifndef QUADRATIC
   if (1) {    

branches/pre/include/ClpModel.hpp

@@ -308 +308 @@
    /// Clp Matrix 
    inline ClpMatrixBase * clpMatrix() const     { return matrix_; }
-   /// Quadratic objective
-   inline CoinPackedMatrix * quadraticObjective() const     { return quadraticObjective_; }
    /// Objective value
    inline double objectiveValue() const {
@@ -488 +486 @@
   /// Row copy if wanted
   ClpMatrixBase * rowCopy_;
-  /// Quadratic objective if any
-  CoinPackedMatrix * quadraticObjective_;
   /// Infeasible/unbounded ray
   double * ray_;

branches/pre/include/ClpSimplex.hpp

@@ -195 +195 @@
       Return code is 0 if nothing interesting, -1 if infeasible both
       ways and +1 if infeasible one way (check values to see which one(s))
+      Solutions are filled in as well - even down, odd up - also
+      status and number of iterations
   */
   int strongBranching(int numberVariables,const int * variables,
                       double * newLower, double * newUpper,
+                      double ** outputSolution,
+                      int * outputStatus, int * outputIterations,
                       bool stopOnFirstInfeasible=true,
                       bool alwaysFinish=false);

branches/pre/include/ClpSimplexDual.hpp

@@ -120 +120 @@
       Return code is 0 if nothing interesting, -1 if infeasible both
       ways and +1 if infeasible one way (check values to see which one(s))
+      Solutions are filled in as well - even down, odd up - also
+      status and number of iterations
   */
   int strongBranching(int numberVariables,const int * variables,
                       double * newLower, double * newUpper,
+                      double ** outputSolution,
+                      int * outputStatus, int * outputIterations,
                       bool stopOnFirstInfeasible=true,
                       bool alwaysFinish=false);

branches/pre/include/ClpSimplexPrimalQuadratic.hpp

@@ -12 +12 @@
 
 class ClpQuadraticInfo;
+class ClpQuadraticObjective;
 
 #include "ClpSimplexPrimal.hpp"
@@ -60 +61 @@
                    ClpQuadraticInfo & info);
   /// Checks complementarity and computes infeasibilities
-  int checkComplementarity (const ClpQuadraticInfo * info,
+  int checkComplementarity (ClpQuadraticInfo * info,
                             CoinIndexedVector * array1, CoinIndexedVector * array2);
   /// Fills in reduced costs
@@ -158 +159 @@
   {return currentSolution_;};
   void setCurrentSolution(const double * solution);
-  /// Original objective
-  inline const double * originalObjective() const
-  {return originalModel_->objective();};
   /// Quadratic sequence or -1 if linear
   inline const int * quadraticSequence() const
@@ -167 +165 @@
   inline const int * backSequence() const
   {return backSequence_;};
-  /// Returns pointer to original model
-  inline const ClpSimplex * originalModel() const
-  { return originalModel_;};
+  /// Returns pointer to original objective
+  inline const ClpQuadraticObjective * originalObjective() const
+  { return originalObjective_;};
+  inline void setOriginalObjective( ClpQuadraticObjective * obj)
+  { originalObjective_ = obj;};
+  /// Quadratic objective
+  CoinPackedMatrix * quadraticObjective() const;
+  /// Linear objective
+  double * linearObjective() const;
   /// Save current In and Sj status
   void saveStatus();
@@ -177 +181 @@
   inline double * djWeight() const
   { return djWeight_;};
+  /// Infeas cost
+  inline double infeasCost() const
+  { return infeasCost_;};
+  inline void setInfeasCost(double value)
+  { infeasCost_ = value;};
   //@}
     
@@ -182 +191 @@
   /**@name Data members */
   //@{
-  /// Model
-  const ClpSimplex * originalModel_;
+  /// Objective
+  const ClpQuadraticObjective * originalObjective_;
   /// Quadratic sequence
   int * quadraticSequence_;
@@ -212 +221 @@
   /// Number of quadratic columns 
   int numberQuadraticColumns_;
+  /// Infeasibility cost
+  double infeasCost_;
   //@}
 };