Changeset 687 for trunk

Timestamp:

Jul 15, 2007 5:28:33 PM (13 years ago)

Author:

forrest

Message:

add threaded to trunk

Location:

trunk/Cbc/src

Files:

• CbcBranchBase.cpp (modified) (2 diffs)
• CbcBranchBase.hpp (modified) (4 diffs)
• CbcBranchDynamic.cpp (modified) (7 diffs)
• CbcBranchDynamic.hpp (modified) (1 diff)
• CbcCompareBase.hpp (modified) (3 diffs)
• CbcCountRowCut.cpp (modified) (4 diffs)
• CbcCutGenerator.cpp (modified) (7 diffs)
• CbcCutGenerator.hpp (modified) (3 diffs)
• CbcFathom.cpp (modified) (2 diffs)
• CbcFathom.hpp (modified) (2 diffs)
• CbcHeuristic.cpp (modified) (5 diffs)
• CbcHeuristic.hpp (modified) (1 diff)
• CbcHeuristicFPump.cpp (modified) (19 diffs)
• CbcHeuristicFPump.hpp (modified) (1 diff)
• CbcHeuristicGreedy.cpp (modified) (1 diff)
• CbcHeuristicLocal.cpp (modified) (1 diff)
• CbcLinked.cpp (modified) (91 diffs)
• CbcLinked.hpp (modified) (18 diffs)
• CbcMessage.cpp (modified) (3 diffs)
• CbcMessage.hpp (modified) (2 diffs)
• CbcModel.cpp (modified) (77 diffs)
• CbcModel.hpp (modified) (7 diffs)
• CbcNode.cpp (modified) (12 diffs)
• CbcSolver.cpp (modified) (55 diffs)
• CbcTree.cpp (modified) (2 diffs)
• Cbc_ampl.cpp (modified) (19 diffs)
• ClpAmplObjective.cpp (added)
• ClpAmplObjective.hpp (added)
• ClpAmplStuff.cpp (added)
• ClpConstraintAmpl.hpp (added)
• Makefile.am (modified) (3 diffs)
• Makefile.in (modified) (5 diffs)

trunk/Cbc/src/CbcBranchBase.cpp

@@ -179 +179 @@
 {
   return NULL;
+}
+/* Pass in information on branch just done and create CbcObjectUpdateData instance.
+   If object does not need data then backward pointer will be NULL.
+   Assumes can get information from solver */
+CbcObjectUpdateData 
+CbcObject::createUpdateInformation(const OsiSolverInterface * solver, 
+                                                        const CbcNode * node,
+                                                        const CbcBranchingObject * branchingObject)
+{
+  return CbcObjectUpdateData();
 }
   
@@ -318 +328 @@
   return *this;
 }
-
-  
+// Default constructor
+CbcObjectUpdateData::CbcObjectUpdateData()
+  : object_(NULL),
+    way_(0),
+    objectNumber_(-1),
+    change_(0.0),
+    status_(0),
+    intDecrease_(0),
+    branchingValue_(0.0)
+{
+}
+
+// Useful constructor
+CbcObjectUpdateData::CbcObjectUpdateData (CbcObject * object,
+                                          int way,
+                                          double change,
+                                          int status,
+                                          int intDecrease,
+                                          double branchingValue)
+  : object_(object),
+    way_(way),
+    objectNumber_(-1),
+    change_(change),
+    status_(status),
+    intDecrease_(intDecrease),
+    branchingValue_(branchingValue)
+{
+}
+
+// Destructor 
+CbcObjectUpdateData::~CbcObjectUpdateData ()
+{
+}
+
+// Copy constructor 
+CbcObjectUpdateData::CbcObjectUpdateData ( const CbcObjectUpdateData & rhs)
+  : object_(rhs.object_),
+    way_(rhs.way_),
+    objectNumber_(rhs.objectNumber_),
+    change_(rhs.change_),
+    status_(rhs.status_),
+    intDecrease_(rhs.intDecrease_),
+    branchingValue_(rhs.branchingValue_)
+{
+}
+
+// Assignment operator 
+CbcObjectUpdateData & 
+CbcObjectUpdateData::operator=( const CbcObjectUpdateData& rhs)
+{
+  if (this!=&rhs) {
+    object_ = rhs.object_;
+    way_ = rhs.way_;
+    objectNumber_ = rhs.objectNumber_;
+    change_ = rhs.change_;
+    status_ = rhs.status_;
+    intDecrease_ = rhs.intDecrease_;
+    branchingValue_ = rhs.branchingValue_;
+  }
+  return *this;
+}
+
+  

trunk/Cbc/src/CbcBranchBase.hpp

@@ -15 +15 @@
 class CbcBranchingObject;
 class OsiChooseVariable;
+class CbcObjectUpdateData;
 
 //#############################################################################
@@ -190 +191 @@
                             double tolerance) const;
 
+  /** Pass in information on branch just done and create CbcObjectUpdateData instance.
+      If object does not need data then backward pointer will be NULL.
+      Assumes can get information from solver */
+  virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver, 
+                                                        const CbcNode * node,
+                                                        const CbcBranchingObject * branchingObject);
+
+  /// Update object by CbcObjectUpdateData
+  virtual void updateInformation(const CbcObjectUpdateData & data) {};
+
   /// Identifier (normally column number in matrix)
   inline int id() const
@@ -323 +334 @@
   {way_=way;};
 
+   /// update model
+  inline void setModel(CbcModel * model)
+  { model_ = model;};
   /// Return model
   inline CbcModel * model() const
@@ -489 +503 @@
 protected:
 };
+/*  This stores data so an object can be updated
+ */
+class CbcObjectUpdateData {
+
+public:
+
+  /// Default Constructor 
+  CbcObjectUpdateData ();
+
+  /// Useful constructor
+  CbcObjectUpdateData (CbcObject * object,
+                       int way,
+                       double change,
+                       int status,
+                       int intDecrease_,
+                       double branchingValue);
+  
+  /// Copy constructor 
+  CbcObjectUpdateData ( const CbcObjectUpdateData &);
+   
+  /// Assignment operator 
+  CbcObjectUpdateData & operator=( const CbcObjectUpdateData& rhs);
+
+  /// Destructor 
+  virtual ~CbcObjectUpdateData ();
+
+  
+public:
+  /// data
+
+  /// Object
+  CbcObject * object_;
+  /// Branch as defined by instance of CbcObject
+  int way_;
+  /// Object number
+  int objectNumber_;
+  /// Change in objective
+  double change_;
+  /// Status 0 Optimal, 1 infeasible, 2 unknown
+  int status_;
+  /// Decrease in number unsatisfied
+  int intDecrease_;
+  /// Branching value
+  double branchingValue_;
+
+};
 
 #endif

trunk/Cbc/src/CbcBranchDynamic.cpp

@@ -267 +267 @@
 {
 }
+// Copy some information i.e. just variable stuff
+void 
+CbcSimpleIntegerDynamicPseudoCost::copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject)
+{
+  downDynamicPseudoCost_=otherObject->downDynamicPseudoCost_;
+  upDynamicPseudoCost_=otherObject->upDynamicPseudoCost_;
+  sumDownCost_ = otherObject->sumDownCost_;
+  sumUpCost_ = otherObject->sumUpCost_;
+  sumDownChange_ = otherObject->sumDownChange_;
+  sumUpChange_ = otherObject->sumUpChange_;
+  sumDownCostSquared_ = otherObject->sumDownCostSquared_;
+  sumUpCostSquared_ = otherObject->sumUpCostSquared_;
+  sumDownDecrease_ = otherObject->sumDownDecrease_;
+  sumUpDecrease_ = otherObject->sumUpDecrease_;
+  lastDownCost_ = otherObject->lastDownCost_;
+  lastUpCost_ = otherObject->lastUpCost_;
+  lastDownDecrease_ = otherObject->lastDownDecrease_;
+  lastUpDecrease_ = otherObject->lastUpDecrease_;
+  numberTimesDown_ = otherObject->numberTimesDown_;
+  numberTimesUp_ = otherObject->numberTimesUp_;
+  numberTimesDownInfeasible_ = otherObject->numberTimesDownInfeasible_;
+  numberTimesUpInfeasible_ = otherObject->numberTimesUpInfeasible_;
+  numberTimesDownLocalFixed_ = otherObject->numberTimesDownLocalFixed_;
+  numberTimesUpLocalFixed_ = otherObject->numberTimesUpLocalFixed_;
+  numberTimesDownTotalFixed_ = otherObject->numberTimesDownTotalFixed_;
+  numberTimesUpTotalFixed_ = otherObject->numberTimesUpTotalFixed_;
+  numberTimesProbingTotal_ = otherObject->numberTimesProbingTotal_;
+}
 // Creates a branching object
 CbcBranchingObject * 
@@ -646 +674 @@
   double downCost = CoinMax((value-below)*downDynamicPseudoCost_,0.0);
   return downCost;
+}
+/* Pass in information on branch just done and create CbcObjectUpdateData instance.
+   If object does not need data then backward pointer will be NULL.
+   Assumes can get information from solver */
+CbcObjectUpdateData 
+CbcSimpleIntegerDynamicPseudoCost::createUpdateInformation(const OsiSolverInterface * solver, 
+                                                           const CbcNode * node,
+                                                           const CbcBranchingObject * branchingObject)
+{
+  double originalValue=node->objectiveValue();
+  int originalUnsatisfied = node->numberUnsatisfied();
+  double objectiveValue = solver->getObjValue()*solver->getObjSense();
+  int unsatisfied=0;
+  int i;
+  //might be base model - doesn't matter
+  int numberIntegers = model_->numberIntegers();;
+  const double * solution = solver->getColSolution();
+  //const double * lower = solver->getColLower();
+  //const double * upper = solver->getColUpper();
+  double change = CoinMax(0.0,objectiveValue-originalValue);
+  int iStatus;
+  if (solver->isProvenOptimal())
+    iStatus=0; // optimal
+  else if (solver->isIterationLimitReached()
+           &&!solver->isDualObjectiveLimitReached())
+    iStatus=2; // unknown 
+  else
+    iStatus=1; // infeasible
+
+  bool feasible = iStatus!=1;
+  if (feasible) {
+    double integerTolerance = 
+      model_->getDblParam(CbcModel::CbcIntegerTolerance);
+    const int * integerVariable = model_->integerVariable();
+    for (i=0;i<numberIntegers;i++) {
+      int j=integerVariable[i];
+      double value = solution[j];
+      double nearest = floor(value+0.5);
+      if (fabs(value-nearest)>integerTolerance) 
+        unsatisfied++;
+    }
+  }
+  int way = branchingObject->way();
+  way = - way; // because after branch so moved on
+  double value = branchingObject->value();
+  return CbcObjectUpdateData (this, way,
+                                  change, iStatus,
+                                  originalUnsatisfied-unsatisfied,value);
+}
+// Update object by CbcObjectUpdateData
+void 
+CbcSimpleIntegerDynamicPseudoCost::updateInformation(const CbcObjectUpdateData & data)
+{
+  bool feasible = data.status_!=1;
+  int way = data.way_;
+  double value = data.branchingValue_;
+  double change = data.change_;
+#define TYPE2 1
+#define TYPERATIO 0.9
+  if (way<0) {
+    // down
+    if (feasible) {
+      //printf("(down change %g value down %g ",change,value-floor(value));
+      incrementNumberTimesDown();
+      addToSumDownChange(1.0e-30+value-floor(value));
+      addToSumDownDecrease(data.intDecrease_);
+#if TYPE2==0
+      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
+      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
+#elif TYPE2==1
+      addToSumDownCost(change);
+      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
+#elif TYPE2==2
+      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
+      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
+#endif
+    } else {
+      //printf("(down infeasible value down %g ",change,value-floor(value));
+      incrementNumberTimesDownInfeasible();
+#if INFEAS==2
+      double distanceToCutoff=0.0;
+      double objectiveValue = model->getCurrentMinimizationObjValue();
+      distanceToCutoff =  model->getCutoff()  - originalValue;
+      if (distanceToCutoff<1.0e20) 
+        change = distanceToCutoff*2.0;
+      else 
+        change = downDynamicPseudoCost()*(value-floor(value))*10.0; 
+      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
+      incrementNumberTimesDown();
+      addToSumDownChange(1.0e-30+value-floor(value));
+      addToSumDownDecrease(data.intDecrease_);
+#if TYPE2==0
+      addToSumDownCost(change/(1.0e-30+(value-floor(value))));
+      setDownDynamicPseudoCost(sumDownCost()/(double) numberTimesDown());
+#elif TYPE2==1
+      addToSumDownCost(change);
+      setDownDynamicPseudoCost(sumDownCost()/sumDownChange());
+#elif TYPE2==2
+      addToSumDownCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(value-floor(value))));
+      setDownDynamicPseudoCost(sumDownCost()*(TYPERATIO/sumDownChange()+(1.0-TYPERATIO)/(double) numberTimesDown()));
+#endif
+#endif
+    }
+  } else {
+    // up
+    if (feasible) {
+      //printf("(up change %g value down %g ",change,ceil(value)-value);
+      incrementNumberTimesUp();
+      addToSumUpChange(1.0e-30+ceil(value)-value);
+      addToSumUpDecrease(data.intDecrease_);
+#if TYPE2==0
+      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
+      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
+#elif TYPE2==1
+      addToSumUpCost(change);
+      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
+#elif TYPE2==2
+      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
+      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
+#endif
+    } else {
+      //printf("(up infeasible value down %g ",change,ceil(value)-value);
+      incrementNumberTimesUpInfeasible();
+#if INFEAS==2
+      double distanceToCutoff=0.0;
+      double objectiveValue = model->getCurrentMinimizationObjValue();
+      distanceToCutoff =  model->getCutoff()  - originalValue;
+      if (distanceToCutoff<1.0e20) 
+        change = distanceToCutoff*2.0;
+      else 
+        change = upDynamicPseudoCost()*(ceil(value)-value)*10.0; 
+      change = CoinMax(1.0e-12*(1.0+fabs(originalValue)),change);
+      incrementNumberTimesUp();
+      addToSumUpChange(1.0e-30+ceil(value)-value);
+      addToSumUpDecrease(data.intDecrease_);
+#if TYPE2==0
+      addToSumUpCost(change/(1.0e-30+(ceil(value)-value)));
+      setUpDynamicPseudoCost(sumUpCost()/(double) numberTimesUp());
+#elif TYPE2==1
+      addToSumUpCost(change);
+      setUpDynamicPseudoCost(sumUpCost()/sumUpChange());
+#elif TYPE2==2
+      addToSumUpCost(change*TYPERATIO+(1.0-TYPERATIO)*change/(1.0e-30+(ceil(value)-value)));
+      setUpDynamicPseudoCost(sumUpCost()*(TYPERATIO/sumUpChange()+(1.0-TYPERATIO)/(double) numberTimesUp()));
+#endif
+#endif
+    }
+  }
+  //print(1,0.5);
 }
 // Pass in probing information
@@ -952 +1129 @@
   int way = object_->way();
   double value = object_->value();
-#define TYPE2 1
-#define TYPERATIO 0.9
+  //#define TYPE2 1
+  //#define TYPERATIO 0.9
   if (way<0) {
     // down
@@ -1043 +1220 @@
     }
   }
-  //object->print();
+  //object->print(1,0.5);
   delete object_;
   object_=NULL;
@@ -1061 +1238 @@
                             double changeDown, int numInfDown)
 {
-  int stateOfSearch = thisOne->model()->stateOfSearch();
+  CbcModel * model = thisOne->model();
+  int stateOfSearch = model->stateOfSearch();
   int betterWay=0;
   double value=0.0;
+  if (!bestObject_) {
+    bestCriterion_=-1.0;
+    bestNumberUp_=INT_MAX;
+    bestNumberDown_=INT_MAX;
+  }
   if (stateOfSearch<=1&&thisOne->model()->currentNode()->depth()>10) {
-#if 0
-    if (!bestObject_) {
-      bestNumberUp_=INT_MAX;
-      bestNumberDown_=INT_MAX;
-    }
+#define TRY_STUFF 1
+#ifdef TRY_STUFF
     // before solution - choose smallest number 
     // could add in depth as well
@@ -1108 +1288 @@
     }
 #else
-    if (!bestObject_) {
-      bestCriterion_=-1.0;
-    }
     // got a solution
     double minValue = CoinMin(changeDown,changeUp);
@@ -1124 +1301 @@
 #endif
   } else {
-    if (!bestObject_) {
-      bestCriterion_=-1.0;
-    }
+#if TRY_STUFF > 1
+    // Get current number of infeasibilities, cutoff and current objective
+    CbcNode * node = model->currentNode();
+    int numberUnsatisfied = node->numberUnsatisfied();
+    double cutoff = model->getCutoff();
+    double objectiveValue = node->objectiveValue();
+#endif
     // got a solution
     double minValue = CoinMin(changeDown,changeUp);
     double maxValue = CoinMax(changeDown,changeUp);
     // Reduce
+#ifdef TRY_STUFF
+    //maxValue = CoinMin(maxValue,minValue*4.0);
+#else
     maxValue = CoinMin(maxValue,minValue*2.0);
+#endif
     value = WEIGHT_AFTER*minValue + (1.0-WEIGHT_AFTER)*maxValue;
-    if (value>bestCriterion_+1.0e-8) {
+    double useValue = value;
+    double useBest = bestCriterion_;
+#if TRY_STUFF>1
+    int thisNumber = CoinMin(numInfUp,numInfDown);
+    int bestNumber = CoinMin(bestNumberUp_,bestNumberDown_);
+    double distance = cutoff-objectiveValue;
+    assert (distance>=0.0);
+    if (useValue+0.1*distance>useBest&&useValue*1.1>useBest&&
+        useBest+0.1*distance>useValue&&useBest*1.1>useValue) {
+      // not much in it - look at unsatisfied
+      if (thisNumber<numberUnsatisfied||bestNumber<numberUnsatisfied) {
+        double perInteger = distance/ ((double) numberUnsatisfied);
+        useValue += thisNumber*perInteger;
+        useBest += bestNumber*perInteger;
+      }
+    }
+#endif
+    if (useValue>useBest+1.0e-8) {
       if (changeUp<=changeDown) {
         betterWay=1;

trunk/Cbc/src/CbcBranchDynamic.hpp

@@ -64 +64 @@
   virtual CbcBranchingObject * createBranch(OsiSolverInterface * solver,
                                             const OsiBranchingInformation * info, int way) ;
+
+  /** Pass in information on branch just done and create CbcObjectUpdateData instance.
+      If object does not need data then backward pointer will be NULL.
+      Assumes can get information from solver */
+  virtual CbcObjectUpdateData createUpdateInformation(const OsiSolverInterface * solver, 
+                                                        const CbcNode * node,
+                                                        const CbcBranchingObject * branchingObject);
+  /// Update object by CbcObjectUpdateData
+  virtual void updateInformation(const CbcObjectUpdateData & data) ;
+  /// Copy some information i.e. just variable stuff
+  void copySome(CbcSimpleIntegerDynamicPseudoCost * otherObject);
 
   /** Create an OsiSolverBranch object

trunk/Cbc/src/CbcCompareBase.hpp

@@ -18 +18 @@
 */
 #include "CbcNode.hpp"
+#include "CbcConfig.h"
 
 class CbcModel;
@@ -76 +77 @@
   {
     assert (x);
+    assert (y);
+#ifndef CBC_THREAD
     CbcNodeInfo * infoX = x->nodeInfo();
     assert (infoX);
     int nodeNumberX = infoX->nodeNumber();
-    assert (y);
     CbcNodeInfo * infoY = y->nodeInfo();
     assert (infoY);
@@ -85 +87 @@
     assert (nodeNumberX!=nodeNumberY);
     return (nodeNumberX>nodeNumberY);
+#else
+    // doesn't work if threaded
+    assert (x!=y);
+    return (x>y);
+#endif
   };
 protected:

trunk/Cbc/src/CbcCountRowCut.cpp

@@ -41 +41 @@
                                 int whichGenerator)
   : OsiRowCut(rhs),
-  owner_(info),
-  ownerCut_(whichOne),
+    owner_(info),
+    ownerCut_(whichOne),
     numberPointingToThis_(0),
-  whichCutGenerator_(whichGenerator)
+    whichCutGenerator_(whichGenerator)
 {
 #ifdef CHECK_CUT_COUNTS
@@ -58 +58 @@
   // Look at owner and delete
   owner_->deleteCut(ownerCut_);
+  ownerCut_=-1234567;
 }
 // Increment number of references
@@ -63 +64 @@
 CbcCountRowCut::increment(int change)
 {
+  assert(ownerCut_!=-1234567);
   numberPointingToThis_+=change;
 }
@@ -70 +72 @@
 CbcCountRowCut::decrement(int change)
 {
-  assert(numberPointingToThis_>=change);
+  assert(ownerCut_!=-1234567);
+  //assert(numberPointingToThis_>=change);
+  assert(numberPointingToThis_>=0);
+  if(numberPointingToThis_<change) {
+    assert(numberPointingToThis_>0);
+    printf("negative cut count %d - %d\n",numberPointingToThis_, change);
+    change = numberPointingToThis_;
+  }
   numberPointingToThis_-=change;
   return numberPointingToThis_;
 }
+
 // Set information
 void 

trunk/Cbc/src/CbcCutGenerator.cpp

@@ -83 +83 @@
   model_ = rhs.model_;
   generator_=rhs.generator_->clone();
-  generator_->refreshSolver(model_->solver());
+  //generator_->refreshSolver(model_->solver());
   whenCutGenerator_=rhs.whenCutGenerator_;
   whenCutGeneratorInSub_ = rhs.whenCutGeneratorInSub_;
@@ -160 +160 @@
 */
 bool
-CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, CbcNode * node)
+CbcCutGenerator::generateCuts( OsiCuts & cs , bool fullScan, OsiSolverInterface * solver, CbcNode * node)
 {
   int howOften = whenCutGenerator_;
@@ -172 +172 @@
     howOften=1;
   bool returnCode=false;
-  OsiSolverInterface * solver = model_->solver();
+  //OsiSolverInterface * solver = model_->solver();
   int depth;
   if (node)
@@ -210 +210 @@
     if (!generator) {
       // Pass across model information in case it could be useful
-      //OsiSolverInterface * solver = model_->solver();
       //void * saveData = solver->getApplicationData();
       //solver->setApplicationData(model_);
@@ -235 +234 @@
       int numberColumns = solver->getNumCols();
       double primalTolerance = 1.0e-8;
-#if 0
-      int numberChanged=0,ifCut=0;
-      CoinPackedVector lbs;
-      CoinPackedVector ubs;
-      for (j=0;j<numberColumns;j++) {
-        if (solver->isInteger(j)) {
-          if (tightUpper[j]<upper[j]) {
-            numberChanged++;
-            assert (tightUpper[j]==floor(tightUpper[j]+0.5));
-            ubs.insert(j,tightUpper[j]);
-            if (tightUpper[j]<solution[j]-primalTolerance)
-              ifCut=1;
-          }
-          if (tightLower[j]>lower[j]) {
-            numberChanged++;
-            assert (tightLower[j]==floor(tightLower[j]+0.5));
-            lbs.insert(j,tightLower[j]);
-            if (tightLower[j]>solution[j]+primalTolerance)
-              ifCut=1;
-          }
-        } else {
-          if (tightUpper[j]==tightLower[j]&&
-              upper[j]>lower[j]) {
-            // fix
-            //solver->setColLower(j,tightLower[j]);
-            //solver->setColUpper(j,tightUpper[j]);
-            double value = tightUpper[j];
-            numberChanged++;
-            if (value<upper[j])
-              ubs.insert(j,value);
-            if (value>lower[j])
-              lbs.insert(j,value);
-          }
-        }
-      }
-      if (numberChanged) {
-        OsiColCut cc;
-        cc.setUbs(ubs);
-        cc.setLbs(lbs);
-        if (ifCut) {
-          cc.setEffectiveness(100.0);
-        } else {
-          cc.setEffectiveness(1.0e-5);
-        }
-        cs.insert(cc);
-      }
-      // need to resolve if some bounds changed
-      returnCode = !solver->basisIsAvailable();
-      assert (!returnCode);
-#else
       const char * tightenBounds = generator->tightenBounds();
-      for (j=0;j<numberColumns;j++) {
-        if (solver->isInteger(j)) {
-          if (tightUpper[j]<upper[j]) {
-            double nearest = floor(tightUpper[j]+0.5);
-            //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
-            solver->setColUpper(j,nearest);
-            if (nearest<solution[j]-primalTolerance)
-              returnCode=true;
-          }
-          if (tightLower[j]>lower[j]) {
-            double nearest = floor(tightLower[j]+0.5);
-            //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
-            solver->setColLower(j,nearest);
-            if (nearest>solution[j]+primalTolerance)
-              returnCode=true;
-          }
-        } else {
-          if (upper[j]>lower[j]) {
-            if (tightUpper[j]==tightLower[j]) {
-              // fix
-              solver->setColLower(j,tightLower[j]);
-              solver->setColUpper(j,tightUpper[j]);
-              if (tightLower[j]>solution[j]+primalTolerance||
-                  tightUpper[j]<solution[j]-primalTolerance)
-                returnCode=true;
-            } else if (tightenBounds&&tightenBounds[j]) {
-              solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
-              solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
-              if (tightLower[j]>solution[j]+primalTolerance||
-                  tightUpper[j]<solution[j]-primalTolerance)
+      if ((model_->getThreadMode()&2)==0) {
+        for (j=0;j<numberColumns;j++) {
+          if (solver->isInteger(j)) {
+            if (tightUpper[j]<upper[j]) {
+              double nearest = floor(tightUpper[j]+0.5);
+              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
+              solver->setColUpper(j,nearest);
+              if (nearest<solution[j]-primalTolerance)
                 returnCode=true;
             }
+            if (tightLower[j]>lower[j]) {
+              double nearest = floor(tightLower[j]+0.5);
+              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
+              solver->setColLower(j,nearest);
+              if (nearest>solution[j]+primalTolerance)
+                returnCode=true;
+            }
+          } else {
+            if (upper[j]>lower[j]) {
+              if (tightUpper[j]==tightLower[j]) {
+                // fix
+                solver->setColLower(j,tightLower[j]);
+                solver->setColUpper(j,tightUpper[j]);
+                if (tightLower[j]>solution[j]+primalTolerance||
+                    tightUpper[j]<solution[j]-primalTolerance)
+                  returnCode=true;
+              } else if (tightenBounds&&tightenBounds[j]) {
+                solver->setColLower(j,CoinMax(tightLower[j],lower[j]));
+                solver->setColUpper(j,CoinMin(tightUpper[j],upper[j]));
+                if (tightLower[j]>solution[j]+primalTolerance||
+                    tightUpper[j]<solution[j]-primalTolerance)
+                  returnCode=true;
+              }
+            }
           }
-        }
+        }
+      } else {
+        CoinPackedVector lbs;
+        CoinPackedVector ubs;
+        int numberChanged=0;
+        bool ifCut=false;
+        for (j=0;j<numberColumns;j++) {
+          if (solver->isInteger(j)) {
+            if (tightUpper[j]<upper[j]) {
+              double nearest = floor(tightUpper[j]+0.5);
+              //assert (fabs(tightUpper[j]-nearest)<1.0e-5); may be infeasible
+              ubs.insert(j,nearest);
+              numberChanged++;
+              if (nearest<solution[j]-primalTolerance)
+                ifCut=true;
+            }
+            if (tightLower[j]>lower[j]) {
+              double nearest = floor(tightLower[j]+0.5);
+              //assert (fabs(tightLower[j]-nearest)<1.0e-5); may be infeasible
+              lbs.insert(j,nearest);
+              numberChanged++;
+              if (nearest>solution[j]+primalTolerance)
+                ifCut=true;
+            }
+          } else {
+            if (upper[j]>lower[j]) {
+              if (tightUpper[j]==tightLower[j]) {
+                // fix
+                lbs.insert(j,tightLower[j]);
+                ubs.insert(j,tightUpper[j]);
+                if (tightLower[j]>solution[j]+primalTolerance||
+                    tightUpper[j]<solution[j]-primalTolerance)
+                  ifCut=true;
+              } else if (tightenBounds&&tightenBounds[j]) {
+                lbs.insert(j,CoinMax(tightLower[j],lower[j]));
+                ubs.insert(j,CoinMin(tightUpper[j],upper[j]));
+                if (tightLower[j]>solution[j]+primalTolerance||
+                    tightUpper[j]<solution[j]-primalTolerance)
+                  ifCut=true;
+              }
+            }
+          }
+        }
+        if (numberChanged) {
+          OsiColCut cc;
+          cc.setUbs(ubs);
+          cc.setLbs(lbs);
+          if (ifCut) {
+            cc.setEffectiveness(100.0);
+          } else {
+            cc.setEffectiveness(1.0e-5);
+          }
+          cs.insert(cc);
+        }
       }
       //if (!solver->basisIsAvailable()) 
       //returnCode=true;
-#endif
 #if 0
       // Pass across info to pseudocosts
@@ -352 +355 @@
       int k ;
       int nOdd=0;
-      const OsiSolverInterface * solver = model_->solver();
+      //const OsiSolverInterface * solver = model_->solver();
       for (k = numberRowCutsAfter-1;k>=numberRowCutsBefore;k--) {
         OsiRowCut & thisCut = cs.rowCut(k) ;
@@ -382 +385 @@
       for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
         OsiRowCut thisCut = cs.rowCut(k) ;
+        if (thisCut.lb()>thisCut.ub()||
+            thisCut.lb()>1.0e8||
+            thisCut.ub()<-1.0e8)
+          printf("cut from %s has bounds %g and %g!\n",
+                 generatorName_,thisCut.lb(),thisCut.ub());
         if (thisCut.lb()<=thisCut.ub()) {
           /* check size of elements.

trunk/Cbc/src/CbcCutGenerator.hpp

@@ -64 +64 @@
     If node then can find out depth
   */
-  bool generateCuts( OsiCuts &cs, bool fullScan,CbcNode * node); 
+  bool generateCuts( OsiCuts &cs, bool fullScan,OsiSolverInterface * solver,CbcNode * node); 
   //@}
 
@@ -177 +177 @@
   inline double timeInCutGenerator() const
   { return timeInCutGenerator_;};
+  inline void incrementTimeInCutGenerator(double value)
+  { timeInCutGenerator_ += value;};
   /// Get the \c CglCutGenerator corresponding to this \c CbcCutGenerator.
   inline CglCutGenerator * generator() const
@@ -237 +239 @@
   inline void setNumberActiveCutsAtRoot(int value)
   { numberActiveCutsAtRoot_ = value;};
+  /// Set model
+  inline void setModel(CbcModel * model)
+  { model_ = model;};
   //@}
   

trunk/Cbc/src/CbcFathom.cpp

@@ -9 +9 @@
 #include <cfloat>
 
-#include "OsiSolverInterface.hpp"
+#ifdef COIN_HAS_CLP
+#include "OsiClpSolverInterface.hpp"
+#endif
 #include "CbcModel.hpp"
 #include "CbcMessage.hpp"
@@ -45 +47 @@
   model_ = model;
 }
+#ifdef COIN_HAS_CLP
+
+//#############################################################################
+// Constructors, destructors clone and assignment
+//#############################################################################
+
+//-------------------------------------------------------------------
+// Default Constructor 
+//-------------------------------------------------------------------
+CbcOsiSolver::CbcOsiSolver ()
+  : OsiClpSolverInterface()
+{
+  cbcModel_ = NULL;
+}
+//-------------------------------------------------------------------
+// Clone
+//-------------------------------------------------------------------
+OsiSolverInterface * 
+CbcOsiSolver::clone(bool copyData) const
+{
+  assert (copyData);
+  return new CbcOsiSolver(*this);
+}
+
+
+//-------------------------------------------------------------------
+// Copy constructor 
+//-------------------------------------------------------------------
+CbcOsiSolver::CbcOsiSolver (
+                  const CbcOsiSolver & rhs)
+  : OsiClpSolverInterface(rhs)
+{
+  cbcModel_ = rhs.cbcModel_;
+}
+
+//-------------------------------------------------------------------
+// Destructor 
+//-------------------------------------------------------------------
+CbcOsiSolver::~CbcOsiSolver ()
+{
+}
+
+//-------------------------------------------------------------------
+// Assignment operator 
+//-------------------------------------------------------------------
+CbcOsiSolver &
+CbcOsiSolver::operator=(const CbcOsiSolver& rhs)
+{
+  if (this != &rhs) { 
+    OsiClpSolverInterface::operator=(rhs);
+    cbcModel_ = rhs.cbcModel_;
+  }
+  return *this;
+}
+#endif

trunk/Cbc/src/CbcFathom.hpp

@@ -3 +3 @@
 #ifndef CbcFathom_H
 #define CbcFathom_H
+#include "CbcConfig.h"
 
 class CbcModel;
@@ -63 +64 @@
   
 };
+#ifdef COIN_HAS_CLP
+#include "OsiClpSolverInterface.hpp"
 
+//#############################################################################
+
+/**
+   
+This is for codes where solver needs to know about CbcModel
+*/
+
+class CbcOsiSolver : public OsiClpSolverInterface {
+  
+public:
+  
+  /**@name Constructors and destructors */
+  //@{
+  /// Default Constructor
+  CbcOsiSolver ();
+  
+  /// Clone
+  virtual OsiSolverInterface * clone(bool copyData=true) const;
+  
+  /// Copy constructor 
+  CbcOsiSolver (const CbcOsiSolver &);
+  
+  /// Assignment operator 
+  CbcOsiSolver & operator=(const CbcOsiSolver& rhs);
+  
+  /// Destructor 
+  virtual ~CbcOsiSolver ();
+  
+  //@}
+  
+  
+  /**@name Sets and Gets */
+  //@{
+  /// Set Cbc Model
+  inline void setCbcModel(CbcModel * model)
+  { cbcModel_=model;};
+  /// Return Cbc Model
+  inline CbcModel * cbcModel() const
+  { return cbcModel_;};
+  //@}
+  
+  //---------------------------------------------------------------------------
+  
+protected:
+  
+  
+  /**@name Private member data */
+  //@{
+  /// Pointer back to CbcModel
+  CbcModel * cbcModel_;
+  //@}
+};
 #endif
+#endif

trunk/Cbc/src/CbcHeuristic.cpp

@@ -178 +178 @@
       // probably faster to use a basis to get integer solutions
       model.setSpecialOptions(2);
+#ifdef CBC_THREAD
+      if (model_->getNumberThreads()>0&&(model_->getThreadMode()&1)!=0) {
+        // See if at root node
+        bool atRoot = model_->getNodeCount()==0;
+        int passNumber = model_->getCurrentPassNumber();
+        if (atRoot&&passNumber==1)
+          model.setNumberThreads(model_->getNumberThreads());
+      }
+#endif
       model.branchAndBound();
       if (logLevel>1)
@@ -302 +311 @@
       (when()%10==2&&(model_->phase()!=2&&model_->phase()!=3)))
     return 0; // switched off
-
   OsiSolverInterface * solver = model_->solver();
   const double * lower = solver->getColLower();
@@ -315 +323 @@
 
   int numberRows = matrix_.getNumRows();
-
+  assert (numberRows<=solver->getNumRows());
   int numberIntegers = model_->numberIntegers();
   const int * integerVariable = model_->integerVariable();
@@ -322 +330 @@
   double newSolutionValue = direction*solver->getObjValue();
   int returnCode = 0;
-
   // Column copy
   const double * element = matrix_.getElements();
@@ -893 +900 @@
       }
     }
+    // Just in case of some stupidity
+    double objOffset=0.0;
+    solver->getDblParam(OsiObjOffset,objOffset);
+    newSolutionValue = -objOffset;
+    for ( i=0 ; i<numberColumns ; i++ )
+      newSolutionValue += objective[i]*newSolution[i];
+    newSolutionValue *= direction;
+    //printf("new solution value %g %g\n",newSolutionValue,solutionValue);
     if (newSolutionValue<solutionValue) {
       // paranoid check

trunk/Cbc/src/CbcHeuristic.hpp

@@ -79 +79 @@
   inline int numberNodes() const
   { return numberNodes_;};
+  /// Just set model - do not do anything else
+  inline void setModelOnly(CbcModel * model)
+  { model_ = model;};
+  
 
   /// Sets fraction of new(rows+columns)/old(rows+columns) before doing small branch and bound (default 1.0)

trunk/Cbc/src/CbcHeuristicFPump.cpp

@@ -175 +175 @@
                          double * betterSolution)
 {
+#define LEN_PRINT 132
+  char pumpPrint[LEN_PRINT];
+  pumpPrint[0]='\0';
   if (!when()||(when()==1&&model_->phase()!=1))
     return 0; // switched off
@@ -364 +367 @@
       }
       if (maximumTime_>0.0&&CoinCpuTime()>=startTime_+maximumTime_) break;
-      numberPasses++;
       memcpy(newSolution,solution,numberColumns*sizeof(double));
       int flip;
-      returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
-                          roundExpensive_,defaultRounding_,&flip);
+      returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
+                          pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
+      numberPasses++;
       if (returnCode) {
         // SOLUTION IS INTEGER
@@ -382 +385 @@
           newSolutionValue += saveObjective[i]*newSolution[i];
         newSolutionValue *= direction;
-        if (model_->logLevel())
-          printf(" - solution found of %g",newSolutionValue);
+        sprintf(pumpPrint+strlen(pumpPrint)," - solution found of %g",newSolutionValue);
         newLineNeeded=false;
         if (newSolutionValue<solutionValue) {
@@ -416 +418 @@
             solutionFound=true;
             if (general&&saveValue!=newSolutionValue)
-              printf(" - cleaned solution of %g\n",solutionValue);
-            else
-              printf("\n");
+              sprintf(pumpPrint+strlen(pumpPrint)," - cleaned solution of %g",solutionValue);
+            if (pumpPrint[0]!='\0')
+              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+                << pumpPrint
+                <<CoinMessageEol;
+            pumpPrint[0]='\0';
           } else {
-          if (model_->logLevel()) 
-            printf(" - not good enough after small branch and bound\n");
+            sprintf(pumpPrint+strlen(pumpPrint)," - not good enough after mini branch and bound");
+            model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+              << pumpPrint
+              <<CoinMessageEol;
+            pumpPrint[0]='\0';
           }
         } else {
-          if (model_->logLevel()) 
-            printf(" - not good enough\n");
+          sprintf(pumpPrint+strlen(pumpPrint)," - not good enough");
+          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+            << pumpPrint
+            <<CoinMessageEol;
+          pumpPrint[0]='\0';
           newLineNeeded=false;
           returnCode=0;
@@ -448 +459 @@
         if (matched || numberPasses%100 == 0) {
           // perturbation
-          if (model_->logLevel())
-            printf("Perturbation applied");
+          sprintf(pumpPrint+strlen(pumpPrint)," perturbation applied");
           newLineNeeded=true;
           for (i=0;i<numberIntegers;i++) {
@@ -531 +541 @@
           memcpy(newSolution,solution,numberColumns*sizeof(double));
           int flip;
-          returnCode = rounds(solver,newSolution,saveObjective,numberIntegers,integerVariable,
-                              roundExpensive_,defaultRounding_,&flip);
+          returnCode = rounds(solver, newSolution,saveObjective,numberIntegers,integerVariable,
+                              pumpPrint,numberPasses,roundExpensive_,defaultRounding_,&flip);
+          numberPasses++;
           if (returnCode) {
             // solution - but may not be better
@@ -540 +551 @@
               newSolutionValue += saveObjective[i]*newSolution[i];
             newSolutionValue *= direction;
-            if (model_->logLevel())
-              printf(" - intermediate solution found of %g",newSolutionValue);
+            sprintf(pumpPrint+strlen(pumpPrint)," - intermediate solution found of %g",newSolutionValue);
             if (newSolutionValue<solutionValue) {
               memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
@@ -636 +646 @@
           }
         }
-        if (model_->logLevel())
-          printf("\npass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
+        if (pumpPrint[0]!='\0')
+          model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+            << pumpPrint
+            <<CoinMessageEol;
+        pumpPrint[0]='\0';
+        sprintf(pumpPrint+strlen(pumpPrint),"Pass %3d: obj. %10.5f --> ", numberPasses+totalNumberPasses,solver->getObjValue());
         if (closestSolution&&solver->getObjValue()<closestObjectiveValue) {
           int i;
@@ -656 +670 @@
       scaleFactor *= weightFactor_;
     } // END WHILE
-    if (newLineNeeded&&model_->logLevel())
-      printf(" - no solution found\n");
+    if (!solutionFound) {
+      sprintf(pumpPrint+strlen(pumpPrint),"No solution found this major pass");
+      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+        << pumpPrint
+        <<CoinMessageEol;
+      pumpPrint[0]='\0';
+    }
     delete solver;
     for ( j=0;j<NUMBER_OLD;j++) 
@@ -724 +743 @@
       newSolver->initialSolve();
       if (!newSolver->isProvenOptimal()) {
-        newSolver->writeMps("bad.mps");
+        //newSolver->writeMps("bad.mps");
         assert (newSolver->isProvenOptimal());
         break;
       }
-      printf("%d integers at bound fixed and %d continuous",
+      sprintf(pumpPrint+strlen(pumpPrint),"Before mini branch and bound, %d integers at bound fixed and %d continuous",
              nFix,nFixC);
-      if (nFixC2+nFixI==0)
-        printf("\n");
-      else
-        printf(" of which %d were internal integer and %d internal continuous\n",
-             nFixI,nFixC2);
+      if (nFixC2+nFixI!=0)
+        sprintf(pumpPrint+strlen(pumpPrint)," of which %d were internal integer and %d internal continuous",
+                nFixI,nFixC2);
+      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+        << pumpPrint
+        <<CoinMessageEol;
+      pumpPrint[0]='\0';
       double saveValue = newSolutionValue;
       returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
-                                       newSolutionValue,"CbcHeuristicLocalAfterFPump");
+                                       cutoff,"CbcHeuristicLocalAfterFPump");
       if ((returnCode&2)!=0) {
         // could add cut
@@ -743 +764 @@
       }
       if (returnCode) {
-        printf("old sol of %g new of %g\n",saveValue,newSolutionValue);
+        sprintf(pumpPrint+strlen(pumpPrint),"Mini branch and bound improved solution from %g to %g",saveValue,newSolutionValue);
+        model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+          << pumpPrint
+          <<CoinMessageEol;
+        pumpPrint[0]='\0';
         memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
         if (fixContinuous) {
@@ -763 +788 @@
             double value = newSolver->getObjValue()*newSolver->getObjSense();
             if (value<newSolutionValue) {
-              printf("freeing continuous gives a solution of %g\n", value);
+              sprintf(pumpPrint+strlen(pumpPrint),"Freeing continuous variables gives a solution of %g", value);
+              model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+                << pumpPrint
+                <<CoinMessageEol;
+              pumpPrint[0]='\0';
               newSolutionValue=value;
               memcpy(betterSolution,newSolver->getColSolution(),numberColumns*sizeof(double));
             }
           } else {
-            newSolver->writeMps("bad3.mps");
+            //newSolver->writeMps("bad3.mps");
           }
         } 
@@ -786 +815 @@
       double gap = relativeIncrement_*fabs(solutionValue);
       cutoff -= CoinMax(CoinMax(gap,absoluteIncrement_),model_->getCutoffIncrement());
-      printf("round again with cutoff of %g\n",cutoff);
+      sprintf(pumpPrint+strlen(pumpPrint),"Round again with cutoff of %g",cutoff);
+      model_->messageHandler()->message(CBC_FPUMP1,model_->messages())
+        << pumpPrint
+        <<CoinMessageEol;
+      pumpPrint[0]='\0';
       if ((accumulate_&3)<2&&usedColumn)
         memset(usedColumn,0,numberColumns);
-      totalNumberPasses += numberPasses;
+      totalNumberPasses += numberPasses-1;
     } else {
       break;
@@ -817 +850 @@
     newSolver->addRow(numberIntegersOrig,closestSolution,
                       lastSolution,-COIN_DBL_MAX,rhs+10.0);
-    double saveValue = newSolutionValue;
+    //double saveValue = newSolutionValue;
     //newSolver->writeMps("sub");
     int returnCode = smallBranchAndBound(newSolver,numberNodes_,newSolution,newSolutionValue,
@@ -857 +890 @@
                           const double * objective,
                           int numberIntegers, const int * integerVariable,
+                          char * pumpPrint, int & iter,
                           bool roundExpensive, double downValue, int *flip)
 {
@@ -865 +899 @@
   int i;
 
-  static int iter = 0;
   int numberColumns = model_->getNumCols();
   // tmp contains the current obj coefficients 
@@ -905 +938 @@
 
   if (nnv > nn) nnv = nn;
-  if (iter != 0&&model_->logLevel())
-    printf("up = %5d , down = %5d", flip_up, flip_down); fflush(stdout);
+  if (iter != 0)
+    sprintf(pumpPrint+strlen(pumpPrint),"up = %5d , down = %5d", flip_up, flip_down);
   *flip = flip_up + flip_down;
   delete [] tmp;
@@ -913 +946 @@
   const double * columnUpper = solver->getColUpper();
   if (*flip == 0 && iter != 0) {
-    if(model_->logLevel())
-      printf(" -- rand = %4d (%4d) ", nnv, nn);
+    sprintf(pumpPrint+strlen(pumpPrint)," -- rand = %4d (%4d) ", nnv, nn);
      for (i = 0; i < nnv; i++) {
        // was solution[list[i]] = 1. - solution[list[i]]; but does that work for 7>=x>=6
@@ -930 +962 @@
      }
      *flip = nnv;
-  } else if (model_->logLevel()) {
-    printf(" ");
+  } else {
+    //sprintf(pumpPrint+strlen(pumpPrint)," ");
   }
   delete [] list; delete [] val;
-  iter++;
+  //iter++;
     
   const double * rowLower = solver->getRowLower();

trunk/Cbc/src/CbcHeuristicFPump.hpp

@@ -165 +165 @@
   int rounds(OsiSolverInterface * solver,double * solution, const double * objective, 
              int numberIntegers, const int * integerVariable,
+             char * pumpPrint,int & passNumber,
              bool roundExpensive=false,
              double downValue=0.5, int *flip=0);

trunk/Cbc/src/CbcHeuristicGreedy.cpp

@@ -764 +764 @@
       memcpy(betterSolution,newSolution,numberColumns*sizeof(double));
       solutionValue = newSolutionValue;
-      model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())
-        << newSolutionValue
-        << "CbcHeuristicGreedy"<<model_->getCurrentSeconds()
-        <<CoinMessageEol;
       returnCode=1;
     }

trunk/Cbc/src/CbcHeuristicLocal.cpp

@@ -549 +549 @@
         returnCode=1;
         solutionValue = newSolutionValue + bestChange;
-        if (bestChange<-1.0e-1)
-          model_->messageHandler()->message(CBC_HEURISTIC_SOLUTION,model_->messages())
-            << solutionValue
-            << "CbcHeuristicLocal"<<model_->getCurrentSeconds()
-            <<CoinMessageEol;
       } else {
         // bad solution - should not happen so debug if see message

trunk/Cbc/src/CbcLinked.cpp

@@ -3 +3 @@
 #include "CbcConfig.h"
 
-#ifndef COIN_HAS_LINK
-#ifdef COIN_HAS_ASL
-#define COIN_HAS_LINK
-#endif
-#endif
 #include "CoinTime.hpp"
 
@@ -13 +8 @@
 #include "CoinModel.hpp"
 #include "ClpSimplex.hpp"
+#ifdef COIN_HAS_ASL
+#include "ClpAmplObjective.hpp"
+#endif
 // returns jColumn (-2 if linear term, -1 if unknown) and coefficient
 static
@@ -92 +90 @@
   return jColumn;
 }
-#ifdef COIN_HAS_LINK
+#include "ClpQuadraticObjective.hpp"
 #include <cassert>
 #if defined(_MSC_VER)
@@ -105 +103 @@
 #include "ClpPackedMatrix.hpp"
 #include "CoinTime.hpp"
-#include "ClpQuadraticObjective.hpp"
 #include "CbcModel.hpp"
 #include "CbcCutGenerator.hpp"
@@ -180 +177 @@
   //sprintf(temp,"cc%d",iPass);
   //writeMps(temp);
+  //writeMps("tight");
+  //exit(33);
   //printf("wrote cc%d\n",iPass);
   OsiClpSolverInterface::initialSolve();
@@ -185 +184 @@
   if (modelPtr_->status()==0&&(secondaryStatus==2||secondaryStatus==4))
     modelPtr_->cleanup(1);
-  if (!isProvenOptimal())
-    writeMps("yy");
+  //if (!isProvenOptimal())
+  //writeMps("yy");
   if (isProvenOptimal()&&quadraticModel_&&modelPtr_->numberColumns()==quadraticModel_->numberColumns()) {
     // see if qp can get better solution
@@ -233 +232 @@
           int numberGenerators = cbcModel_->numberCutGenerators();
           int iGenerator;
+          cbcModel_->lockThread();
           for (iGenerator=0;iGenerator<numberGenerators;iGenerator++) {
             CbcCutGenerator * generator = cbcModel_->cutGenerator(iGenerator);
@@ -263 +263 @@
             }
           }
+          cbcModel_->unlockThread();
         }
       }
@@ -461 +462 @@
           bestSolution_ = CoinCopyOfArray(modelPtr_->getColSolution(),modelPtr_->getNumCols());
           bestObjectiveValue_ = value;
+          if (maxIts<=10000&&cbcModel_) {
+            OsiSolverLink * solver2 = dynamic_cast<OsiSolverLink *> (cbcModel_->solver());
+            assert (solver2);
+            if (solver2!=this) {
+              // in strong branching - need to store in original solver
+              if (value<solver2->bestObjectiveValue_-1.0e-3) {
+                delete [] solver2->bestSolution_;
+                solver2->bestSolution_ = CoinCopyOfArray(bestSolution_,modelPtr_->getNumCols());
+                solver2->bestObjectiveValue_ = value;
+              }
+            }
+          }
           // If model has stored then add cut (if convex)
           if (cbcModel_&&(specialOptions2_&4)!=0&&quadraticModel_) {
@@ -470 +483 @@
               CglStored * gen2 = dynamic_cast<CglStored *> (gen);
               if (gen2) {
+                cbcModel_->lockThread();
                 // add OA cut
-                double offset;
+                double offset=0.0;
                 int numberColumns = quadraticModel_->numberColumns();
                 double * gradient = new double [numberColumns+1];
@@ -492 +506 @@
                     int yColumn = obj->yColumn();
                     if (xColumn!=yColumn) {
-                      double coefficient = 2.0*obj->coefficient();
+                      double coefficient = /* 2.0* */obj->coefficient();
                       gradient[xColumn] += coefficient*solution[yColumn];
                       gradient[yColumn] += coefficient*solution[xColumn];
@@ -520 +534 @@
                 delete [] gradient;
                 delete [] column;
+                cbcModel_->unlockThread();
                 break;
               }
@@ -634 +649 @@
           delete [] lower2;
           delete [] upper2;
-          if (isProvenOptimal())
-            writeMps("zz");
+          //if (isProvenOptimal())
+          //writeMps("zz");
         }
       }
@@ -652 +667 @@
             CglTemporary * gen2 = dynamic_cast<CglTemporary *> (gen);
             if (gen2) {
+              cbcModel_->lockThread();
               const double * solution = getColSolution();
               // add OA cut
-              double offset;
+              double offset=0.0;
               int numberColumns = quadraticModel_->numberColumns();
               double * gradient = new double [numberColumns+1];
@@ -677 +693 @@
                   int yColumn = obj->yColumn();
                   if (xColumn!=yColumn) {
-                    double coefficient = 2.0*obj->coefficient();
+                    double coefficient = /* 2.0* */obj->coefficient();
                     gradient[xColumn] += coefficient*solution[yColumn];
                     gradient[yColumn] += coefficient*solution[xColumn];
@@ -710 +726 @@
               delete [] gradient;
               delete [] column;
+              cbcModel_->unlockThread();
               break;
             }
@@ -735 +752 @@
             //const double * columnLower = modelPtr_->columnLower();
             //const double * columnUpper = modelPtr_->columnUpper();
+            cbcModel_->lockThread();
             for (int iNon=0;iNon<numberNonLinearRows_;iNon++) {
               int iRow = rowNonLinear_[iNon];
@@ -754 +772 @@
                 int yColumn = obj->yColumn();
                 if (xColumn!=yColumn) {
-                  double coefficient = 2.0*obj->coefficient();
+                  double coefficient = /* 2.0* */obj->coefficient();
                   gradient[xColumn] += coefficient*solution[yColumn];
                   gradient[yColumn] += coefficient*solution[xColumn];
@@ -791 +809 @@
                 gen2->addCut(offset-1.0e-7,COIN_DBL_MAX,n,column,gradient);
             }
+            cbcModel_->unlockThread();
             delete [] gradient;
             delete [] column;
@@ -812 +831 @@
 //-------------------------------------------------------------------
 OsiSolverLink::OsiSolverLink ()
-  : OsiClpSolverInterface()
+  : CbcOsiSolver()
 {
   gutsOfDestructor(true);
@@ -893 +912 @@
 */
 OsiSolverLink::OsiSolverLink ( CoinModel & coinModel)
-  : OsiClpSolverInterface()
+  : CbcOsiSolver()
 {
   gutsOfDestructor(true);
@@ -899 +918 @@
 }
 // need bounds
-static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue)
+static void fakeBounds(OsiSolverInterface * solver,int column,double maximumValue,
+                       CoinModel * model1, CoinModel * model2)
 {
   double lo = solver->getColLower()[column];
-  if (lo<-maximumValue)
+  if (lo<-maximumValue) {
     solver->setColLower(column,-maximumValue);
+    if (model1)
+      model1->setColLower(column,-maximumValue);
+    if (model2)
+      model2->setColLower(column,-maximumValue);
+  }
   double up = solver->getColUpper()[column];
-  if (up>maximumValue)
+  if (up>maximumValue) {
     solver->setColUpper(column,maximumValue);
-}
-void OsiSolverLink::load ( CoinModel & coinModel, bool tightenBounds,int logLevel)
+    if (model1)
+      model1->setColUpper(column,maximumValue);
+    if (model2)
+      model2->setColUpper(column,maximumValue);
+  }
+}
+void OsiSolverLink::load ( CoinModel & coinModelOriginal, bool tightenBounds,int logLevel)
 {
   // first check and set up arrays
-  int numberColumns = coinModel.numberColumns();
-  int numberRows = coinModel.numberRows();
+  int numberColumns = coinModelOriginal.numberColumns();
+  int numberRows = coinModelOriginal.numberRows();
   // List of nonlinear entries
   int * which = new int[numberColumns];
@@ -920 +950 @@
   int numberErrors=0;
   for (iColumn=0;iColumn<numberColumns;iColumn++) {
-    CoinModelLink triple=coinModel.firstInColumn(iColumn);
+    CoinModelLink triple=coinModelOriginal.firstInColumn(iColumn);
     bool linear=true;
     int n=0;
     // See if quadratic objective
-    const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
+    const char * expr = coinModelOriginal.getColumnObjectiveAsString(iColumn);
     if (strcmp(expr,"Numeric")) {
       linear=false;
@@ -930 +960 @@
     while (triple.row()>=0) {
       int iRow = triple.row();
-      const char * expr = coinModel.getElementAsString(iRow,iColumn);
+      const char * expr = coinModelOriginal.getElementAsString(iRow,iColumn);
       if (strcmp(expr,"Numeric")) {
         linear=false;
       }
-      triple=coinModel.next(triple);
+      triple=coinModelOriginal.next(triple);
       n++;
     }
@@ -944 +974 @@
   if (!numberVariables_) {
     delete [] which;
-    coinModel_ = coinModel;
+    coinModel_ = coinModelOriginal;
     int nInt=0;
     for (iColumn=0;iColumn<numberColumns;iColumn++) {   
@@ -951 +981 @@
     }
     printf("There are %d integers\n",nInt);
-    loadFromCoinModel(coinModel,true);
+    loadFromCoinModel(coinModelOriginal,true);
     OsiObject ** objects = new OsiObject * [nInt];
     nInt=0;
@@ -968 +998 @@
     return;
   } else {
-    coinModel_ = coinModel;
+    coinModel_ = coinModelOriginal;
     // arrays for tightening bounds
     int * freeRow = new int [numberRows];
@@ -1004 +1034 @@
           ifFirst=false;
         }
-        coinModel.setObjective(iColumn,linearTerm);
+        coinModelOriginal.setObjective(iColumn,linearTerm);
       }
     }
@@ -1039 +1069 @@
             ifFirst=false;
           }
-          coinModel.setElement(iRow,iColumn,linearTerm);
+          coinModelOriginal.setElement(iRow,iColumn,linearTerm);
         }
         triple=coinModel_.next(triple);
@@ -1053 +1083 @@
     }
     printf("There are %d bilinear and %d integers\n",nBi,nInt);
-    loadFromCoinModel(coinModel,true);
+    loadFromCoinModel(coinModelOriginal,true);
+    CoinModel coinModel = coinModelOriginal;
     if (tightenBounds&&numberColumns<100) {
       // first fake bounds
       for (iColumn=0;iColumn<numberColumns;iColumn++) {
         if (tryColumn[iColumn]) {
-          fakeBounds(this,iColumn,defaultBound_);
+          fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
         }
       }
@@ -1091 +1122 @@
             columnLower[iColumn]=value;
             coinModel_.setColumnLower(iColumn,value);
+            coinModel.setColumnLower(iColumn,value);
           }
           objective[iColumn]=-1.0;
@@ -1102 +1134 @@
             columnUpper[iColumn]=value;
             coinModel_.setColumnUpper(iColumn,value);
+            coinModel.setColumnUpper(iColumn,value);
           }
           objective[iColumn]=0.0;
@@ -1127 +1160 @@
       // add in objective as constraint
       objectiveVariable_= numberColumns;
-      objectiveRow_ = modelPtr_->numberRows();
-      addCol(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
+      objectiveRow_ = coinModel.numberRows();
+      coinModel.addColumn(0,NULL,NULL,-COIN_DBL_MAX,COIN_DBL_MAX,1.0);
       int * column = new int[numberColumns+1];
       double * element = new double[numberColumns+1];
-      double * objective = modelPtr_->objective();
+      double * objective = coinModel.objectiveArray();
       int n=0;
-      int starts[2]={0,0};
       for (int i=0;i<numberColumns;i++) {
         if (objective[i]) {
@@ -1143 +1175 @@
       column[n]=objectiveVariable_;
       element[n++]=-1.0;
-      starts[1]=n;
-      double offset = - modelPtr_->objectiveOffset();
+      double offset = - coinModel.objectiveOffset();
       assert (!offset); // get sign right if happens
-      modelPtr_->setObjectiveOffset(0.0);
+      coinModel.setObjectiveOffset(0.0);
       double lowerBound = -COIN_DBL_MAX;
-      addRows(1,starts,column,element,&lowerBound,&offset);
+      coinModel.addRow(n,column,element,lowerBound,offset);
       delete [] column;
       delete [] element;
@@ -1159 +1190 @@
     double * sort = new double [nBi];
     nBi=nInt;
+    const OsiObject ** justBi = const_cast<const OsiObject **> (objects+nInt);
     for (iColumn=0;iColumn<numberColumns;iColumn++) {
       if (quadraticObjective)
@@ -1166 +1198 @@
       if (strcmp(expr,"Numeric")) {
         // need bounds
-        fakeBounds(this,iColumn,defaultBound_);
+        fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
         // value*x*y
         char temp[20000];
@@ -1179 +1211 @@
             if (quadraticObjective) {
               columnQuadratic[numberQuadratic]=jColumn;
-              elementQuadratic[numberQuadratic++]=2.0*value; // convention
+              if (jColumn==iColumn)
+                elementQuadratic[numberQuadratic++]=2.0*value; // convention
+              else
+                elementQuadratic[numberQuadratic++]=1.0*value; // convention
             }
             // need bounds
-            fakeBounds(this,jColumn,defaultBound_);
+            fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
             double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
             if (meshI)
@@ -1209 +1244 @@
               meshJ=0.0;
             }
-            OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
-                                                   nBi-nInt,(const OsiObject **) (objects+nInt));
+            OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,objectiveRow_,value,meshI,meshJ,
+                                                   nBi-nInt,justBi);
             newObj->setPriority(biLinearPriority_);
             objects[nBi++] = newObj;
@@ -1240 +1275 @@
         if (strcmp("Numeric",el)) {
           // need bounds
-          fakeBounds(this,iColumn,defaultBound_);
+          fakeBounds(this,iColumn,defaultBound_,&coinModel,&coinModel_);
           // value*x*y
           char temp[20000];
@@ -1252 +1287 @@
             if (jColumn>=0) {
               // need bounds
-              fakeBounds(this,jColumn,defaultBound_);
+              fakeBounds(this,jColumn,defaultBound_,&coinModel,&coinModel_);
               double meshI = coinModel_.isInteger(iColumn) ? 1.0 : 0.0;
               if (meshI)
@@ -1279 +1314 @@
                 meshJ=0.0;
               }
-              OsiBiLinear * newObj = new OsiBiLinear(this,iColumn,jColumn,iRow,value,meshI,meshJ,
-                                                     nBi-nInt,(const OsiObject **) (objects+nInt));
+              OsiBiLinear * newObj = new OsiBiLinear(&coinModel,iColumn,jColumn,iRow,value,meshI,meshJ,
+                                                     nBi-nInt,justBi);
               newObj->setPriority(biLinearPriority_);
               objects[nBi++] = newObj;
@@ -1308 +1343 @@
              stats[0],stats[1],stats[2],nBi-nInt-nDiff);
     }
+    // reload with all bilinear stuff
+    loadFromCoinModel(coinModel,true);
+    //exit(77);
     nInt=0;
     for (iColumn=0;iColumn<numberColumns;iColumn++) {   
@@ -1458 +1496 @@
     }
   }
+  delete [] which;
+  if ((specialOptions2_&16)!=0)
+    addTighterConstraints();
+}
+// Add reformulated bilinear constraints
+void 
+OsiSolverLink::addTighterConstraints()
+{
+  // This is first attempt - for now get working on trimloss
+  int numberW=0;
+  int * xW = new int[numberObjects_];
+  int * yW = new int[numberObjects_];
+  // Points to firstlambda
+  int * wW = new int[numberObjects_];
+  // Coefficient
+  double * alphaW = new double[numberObjects_];
+  // Objects
+  OsiBiLinear ** objW = new OsiBiLinear * [numberObjects_];
+  int numberColumns = getNumCols();
+  int firstLambda=numberColumns;
+  // set up list (better to rethink and do properly as column ordered)
+  int * list = new int[numberColumns];
+  memset(list,0,numberColumns*sizeof(int));
+  int i;
+  for ( i =0;i<numberObjects_;i++) {
+    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
+    if (obj) {
+      //obj->setBranchingStrategy(4); // ***** temp
+      objW[numberW]=obj;
+      xW[numberW]=obj->xColumn();
+      yW[numberW]=obj->yColumn();
+      list[xW[numberW]]=1;
+      list[yW[numberW]]=1;
+      wW[numberW]=obj->firstLambda();
+      firstLambda = CoinMin(firstLambda,obj->firstLambda());
+      alphaW[numberW]=obj->coefficient();
+      //assert (alphaW[numberW]==1.0); // fix when occurs
+      numberW++;
+    }
+  }
+  int nList = 0;
+  for (i=0;i<numberColumns;i++) {
+    if (list[i])
+      list[nList++]=i;
+  }
+  // set up mark array
+  char * mark = new char [firstLambda*firstLambda];
+  memset(mark,0,firstLambda*firstLambda);
+  for (i=0;i<numberW;i++) {
+    int x = xW[i];
+    int y = yW[i];
+    mark[x*firstLambda+y]=1;
+    mark[y*firstLambda+x]=1;
+  }
+  int numberRows2 = originalRowCopy_->getNumRows();
+  int * addColumn = new int [numberColumns];
+  double * addElement = new double [numberColumns];
+  int * addW = new int [numberColumns];
+  assert (objectiveRow_<0); // fix when occurs
+  for (int iRow=0;iRow<numberRows2;iRow++) {
+    for (int iList=0;iList<nList;iList++) {
+      int kColumn = list[iList];
+      const double * columnLower = getColLower();
+      //const double * columnUpper = getColUpper();
+      const double * rowLower = getRowLower();
+      const double * rowUpper = getRowUpper();
+      const CoinPackedMatrix * rowCopy = getMatrixByRow();
+      const double * element = rowCopy->getElements();
+      const int * column = rowCopy->getIndices();
+      const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
+      const int * rowLength = rowCopy->getVectorLengths();
+      CoinBigIndex j;
+      int numberElements = rowLength[iRow];
+      int n=0;
+      for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
+        int iColumn = column[j];
+        if (iColumn>=firstLambda) {
+          // no good
+          n=-1;
+          break;
+        }
+        if (mark[iColumn*firstLambda+kColumn])
+          n++;
+      }
+      if (n==numberElements) {
+        printf("can add row %d\n",iRow);
+        assert (columnLower[kColumn]>=0); // might be able to fix
+        n=0;
+        for (j=rowStart[iRow];j<rowStart[iRow]+numberElements;j++) {
+          int xColumn=kColumn;
+          int yColumn = column[j];
+          int k;
+          for (k=0;k<numberW;k++) {
+            if ((xW[k]==yColumn&&yW[k]==xColumn)||
+                (yW[k]==yColumn&&xW[k]==xColumn))
+              break;
+          }
+          assert (k<numberW);
+          if (xW[k]!=xColumn) {
+            int temp=xColumn;
+            xColumn=yColumn;
+            yColumn=temp;
+          }
+          addW[n/4]=k;
+          int start = wW[k];
+          double value = element[j];
+          for (int kk=0;kk<4;kk++) {
+            // Dummy value
+            addElement[n]= value;
+            addColumn[n++]=start+kk;
+          }
+        }
+        addColumn[n++] = kColumn;
+        double lo = rowLower[iRow];
+        double up = rowUpper[iRow];
+        if (lo>-1.0e20) {
+          // and tell object
+          for (j=0;j<n-1;j+=4) {
+            int iObject = addW[j/4];
+            objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
+          }
+          addElement[n-1]=-lo;
+          if (lo==up)
+            addRow(n,addColumn,addElement,0.0,0.0);
+          else
+            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
+          matrix_->appendRow(n,addColumn,addElement);
+        }
+        if (up<1.0e20&&up>lo) {
+          // and tell object
+          for (j=0;j<n-1;j+=4) {
+            int iObject = addW[j/4];
+            objW[iObject]->addExtraRow(matrix_->getNumRows(),addElement[j]);
+          }
+          addElement[n-1]=-up;
+          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
+          matrix_->appendRow(n,addColumn,addElement);
+        }
+      }
+    }
+  }
 #if 0
-  //int fake[10]={3,4,2,5,5,3,1,11,2,0};
-  int fake[10]={11,5,5,4,3,3,2,2,1,0};
-  for (int kk=0;kk<10;kk++) {
-    setColUpper(kk,fake[kk]);
-    setColLower(kk,fake[kk]);
+  // possibly do bounds
+  for (int iColumn=0;iColumn<firstLambda;iColumn++) {
+    for (int iList=0;iList<nList;iList++) {
+      int kColumn = list[iList];
+      const double * columnLower = getColLower();
+      const double * columnUpper = getColUpper();
+      if (mark[iColumn*firstLambda+kColumn]) {
+        printf("can add column %d\n",iColumn);
+        assert (columnLower[kColumn]>=0); // might be able to fix
+        int xColumn=kColumn;
+        int yColumn = iColumn;
+        int k;
+        for (k=0;k<numberW;k++) {
+          if ((xW[k]==yColumn&&yW[k]==xColumn)||
+              (yW[k]==yColumn&&xW[k]==xColumn))
+            break;
+        }
+        assert (k<numberW);
+        if (xW[k]!=xColumn) {
+          int temp=xColumn;
+          xColumn=yColumn;
+          yColumn=temp;
+        }
+        int start = wW[k];
+        int n=0;
+        for (int kk=0;kk<4;kk++) {
+          // Dummy value
+          addElement[n]= 1.0e-19;
+          addColumn[n++]=start+kk;
+        }
+        // Tell object about this
+        objW[k]->addExtraRow(matrix_->getNumRows(),1.0);
+        addColumn[n++] = kColumn;
+        double lo = columnLower[iColumn];
+        double up = columnUpper[iColumn];
+        if (lo>-1.0e20) {
+          addElement[n-1]=-lo;
+          if (lo==up)
+            addRow(n,addColumn,addElement,0.0,0.0);
+          else
+            addRow(n,addColumn,addElement,0.0,COIN_DBL_MAX);
+          matrix_->appendRow(n,addColumn,addElement);
+        }
+        if (up<1.0e20&&up>lo) {
+          addElement[n-1]=-up;
+          addRow(n,addColumn,addElement,-COIN_DBL_MAX,0.0);
+          matrix_->appendRow(n,addColumn,addElement);
+        }
+      }
+    }
   }
 #endif
-  delete [] which;
+  delete [] xW;
+  delete [] yW;
+  delete [] wW;
+  delete [] alphaW;
+  delete [] addColumn;
+  delete [] addElement;
+  delete [] addW;
+  delete [] mark;
+  delete [] list;
+  delete [] objW;
 }
 // Set all biLinear priorities on x-x variables
@@ -1491 +1724 @@
         objNew->setYOtherSatisfied(oldSatisfied);
         objNew->setYMeshSize(meshSize);
-        objNew->setXYSatisfied(0.5*meshSize);
+        objNew->setXYSatisfied(0.25*meshSize);
         objNew->setPriority(value);
         objNew->setBranchingStrategy(8);
@@ -1502 +1735 @@
   delete [] newObject;
 }
+/* Set options and priority on all or some biLinear variables
+   1 - on I-I
+   2 - on I-x
+   4 - on x-x
+      or combinations.
+      -1 means leave (for priority value and strategy value)
+*/
+void 
+OsiSolverLink::setBranchingStrategyOnVariables(int strategyValue, int priorityValue,
+                                               int mode)
+{
+  int i;
+  for ( i =0;i<numberObjects_;i++) {
+    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (object_[i]);
+    if (obj) {
+      bool change=false;
+      if (obj->xMeshSize()<1.0&&obj->yMeshSize()<1.0&&(mode&4)!=0)
+        change=true;
+      else if (((obj->xMeshSize()==1.0&&obj->yMeshSize()<1.0)||
+                (obj->xMeshSize()<1.0&&obj->yMeshSize()==1.0))&&(mode&2)!=0)
+        change=true;
+      else if (obj->xMeshSize()==1.0&&obj->yMeshSize()==1.0&&(mode&1)!=0)
+        change=true;
+      else if (obj->xMeshSize()>1.0||obj->yMeshSize()>1.0)
+        abort();
+      if (change) {
+        if (strategyValue>=0)
+          obj->setBranchingStrategy(strategyValue);
+        if (priorityValue>=0)
+          obj->setPriority(priorityValue);
+      }
+    }
+  }
+}
+
 // Say convex (should work it out)
 void 
@@ -1712 +1980 @@
   }
   model.initialSolve();
-  model.writeMps("bad.mps");
+  //model.writeMps("bad.mps");
   // redo number of columns
   numberColumns = model.numberColumns();
@@ -2079 +2347 @@
       char temp[20];
       sprintf(temp,"pass%d.mps",iPass);
-      model.writeMps(temp);
+      //model.writeMps(temp);
 #ifdef CLP_DEBUG
       if (model.status()) {
@@ -2138 +2406 @@
       }
       model.primal(1);
-      model.writeMps("xx.mps");
+      //model.writeMps("xx.mps");
       iPass--;
       goodMove=-1;
@@ -2270 +2538 @@
     cbcModel->cutGenerator(6)->setTiming(true);
     // For now - switch off most heuristics (because CglPreProcess is bad with QP)
-#if 0    
+#if 1    
     CbcHeuristicFPump heuristicFPump(*cbcModel);
     heuristicFPump.setWhen(13);
@@ -2291 +2559 @@
     
     CbcRounding rounding(*cbcModel);
+    rounding.setHeuristicName("rounding");
     cbcModel->addHeuristic(&rounding);
     
@@ -2329 +2598 @@
     // if convex
     if ((specialOptions2()&4)!=0) {
+      if (cbcModel_)
+        cbcModel_->lockThread();
       // add OA cut
       double offset;
@@ -2350 +2621 @@
       delete [] gradient;
       delete [] column;
+      if (cbcModel_)
+        cbcModel_->unlockThread();
     }
     delete qp;
@@ -2369 +2642 @@
   double * solution = CoinCopyOfArray(temp->primalColumnSolution(),numberColumns);
   delete temp;
+  if (mode==0) {
+    return solution;
+  } else if (mode==2) {
+    const double * lower = getColLower();
+    const double * upper = getColUpper();
+    for (int iObject =0;iObject<numberObjects_;iObject++) {
+      OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
+      if (obj&&(obj->priority()<biLinearPriority_||biLinearPriority_<=0)) {
+        int iColumn = obj->columnNumber();
+        double value = solution[iColumn];
+        value = floor(value+0.5);
+        if (fabs(value-solution[iColumn])>0.01) {
+          setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
+          setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,1.0)));
+        } else {
+          // could fix to integer
+          setColLower(iColumn,CoinMax(lower[iColumn],value-CoinMax(defaultBound_,0.0)));
+          setColUpper(iColumn,CoinMin(upper[iColumn],value+CoinMax(defaultBound_,0.0)));
+        }
+      }
+    }
+    return solution;
+  }
   OsiClpSolverInterface newSolver;
   if (mode==1) {
@@ -2376 +2672 @@
     tempModel = coinModel_;
     // solve modified problem
+    char * mark = new char[numberColumns];
+    memset(mark,0,numberColumns);
     for (int iObject =0;iObject<numberObjects_;iObject++) {
       OsiSimpleInteger * obj = dynamic_cast<OsiSimpleInteger *> (object_[iObject]);
@@ -2383 +2681 @@
         value = ceil(value-1.0e-7);
         tempModel.associateElement(coinModel_.columnName(iColumn),value);
-      }
-    }
+        mark[iColumn]=1;
+      }
+      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
+      if (objB) {
+        // if one or both continuous then fix one
+        if (objB->xMeshSize()<1.0) {
+          int xColumn = objB->xColumn();
+          double value = solution[xColumn];
+          tempModel.associateElement(coinModel_.columnName(xColumn),value);
+          mark[xColumn]=1;
+        } else if (objB->yMeshSize()<1.0) {
+          int yColumn = objB->yColumn();
+          double value = solution[yColumn];
+          tempModel.associateElement(coinModel_.columnName(yColumn),value);
+          mark[yColumn]=1;
+        }
+      }
+    }
+    CoinModel * reOrdered = tempModel.reorder(mark);
+    assert (reOrdered);
+    tempModel=*reOrdered;
+    delete reOrdered;
+    delete [] mark;
     newSolver.loadFromCoinModel(tempModel,true);
     for (int iObject =0;iObject<numberObjects_;iObject++) {
@@ -2394 +2713 @@
         newSolver.setColLower(iColumn,value);
         newSolver.setColUpper(iColumn,value);
+      }
+      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (object_[iObject]);
+      if (objB) {
+        // if one or both continuous then fix one
+        if (objB->xMeshSize()<1.0) {
+          int xColumn = objB->xColumn();
+          double value = solution[xColumn];
+          newSolver.setColLower(xColumn,value);
+          newSolver.setColUpper(xColumn,value);
+        } else if (objB->yMeshSize()<1.0) {
+          int yColumn = objB->yColumn();
+          double value = solution[yColumn];
+          newSolver.setColLower(yColumn,value);
+          newSolver.setColUpper(yColumn,value);
+        }
       }
     }
@@ -2491 +2825 @@
     clpModel->setLogLevel(saveLogLevel);
     returnCode=-1; // infeasible//std::cout<<"Problem is infeasible - tightenPrimalBounds!"<<std::endl;
-    clpModel->writeMps("infeas2.mps");
+    //clpModel->writeMps("infeas2.mps");
   } else {
     clpModel->setLogLevel(saveLogLevel);
@@ -2710 +3044 @@
 {
   assert (copyData);
-  return new OsiSolverLink(*this);
+  OsiSolverLink * newModel = new OsiSolverLink(*this);
+  return newModel;
 }
 
@@ -2719 +3054 @@
 OsiSolverLink::OsiSolverLink (
                   const OsiSolverLink & rhs)
-  : OsiClpSolverInterface(rhs)
+  : CbcOsiSolver(rhs)
 {
   gutsOfDestructor(true);
@@ -2743 +3078 @@
   if (this != &rhs) { 
     gutsOfDestructor();
-    OsiClpSolverInterface::operator=(rhs);
+    CbcOsiSolver::operator=(rhs);
     gutsOfCopy(rhs);
   }
@@ -2771 +3106 @@
   convex_ = NULL;
   whichNonLinear_ = NULL;
-  cbcModel_ = NULL;
   info_ = NULL;
   fixVariables_=NULL;
@@ -2801 +3135 @@
   biLinearPriority_ = rhs.biLinearPriority_;
   numberFix_ = rhs.numberFix_;
-  cbcModel_ = rhs.cbcModel_;
   if (numberVariables_) { 
     if (rhs.matrix_)
@@ -3094 +3427 @@
       }
     }
-    newSolver.writeMps("xx");
+    //newSolver.writeMps("xx");
     CbcModel model(newSolver);
     // Now do requested saves and modifications
@@ -3318 +3651 @@
           setRowBounds(i,-COIN_DBL_MAX,COIN_DBL_MAX);
         initialSolve();
-        if (!isProvenOptimal())
-          getModelPtr()->writeMps("bad.mps");
+        //if (!isProvenOptimal())
+        //getModelPtr()->writeMps("bad.mps");
         if (isProvenOptimal()) {
           delete [] bestSolution_;
@@ -4302 +4635 @@
     xyRow_(-1),
     convexity_(-1),
+    numberExtraRows_(0),
+    multiplier_(NULL),
+    extraRow_(NULL),
     chosen_(-1)
 {
@@ -4330 +4666 @@
     xyRow_(xyRow),
     convexity_(-1),
+    numberExtraRows_(0),
+    multiplier_(NULL),
+    extraRow_(NULL),
     chosen_(-1)
 {
@@ -4517 +4856 @@
   }
 }
+// Useful constructor
+OsiBiLinear::OsiBiLinear (CoinModel * coinModel, int xColumn,
+                          int yColumn, int xyRow, double coefficient,
+                          double xMesh, double yMesh,
+                          int numberExistingObjects,const OsiObject ** objects )
+  : OsiObject2(),
+    coefficient_(coefficient),
+    xMeshSize_(xMesh),
+    yMeshSize_(yMesh),
+    xSatisfied_(1.0e-6),
+    ySatisfied_(1.0e-6),
+    xOtherSatisfied_(0.0),
+    yOtherSatisfied_(0.0),
+    xySatisfied_(1.0e-6),
+    xyBranchValue_(0.0),
+    xColumn_(xColumn),
+    yColumn_(yColumn),
+    firstLambda_(-1),
+    branchingStrategy_(0),
+    boundType_(0),
+    xRow_(-1),
+    yRow_(-1),
+    xyRow_(xyRow),
+    convexity_(-1),
+    numberExtraRows_(0),
+    multiplier_(NULL),
+    extraRow_(NULL),
+    chosen_(-1)
+{
+  double columnLower[4];
+  double columnUpper[4];
+  double objective[4];
+  double rowLower[3];
+  double rowUpper[3];
+  CoinBigIndex starts[5];
+  int index[16];
+  double element[16];
+  int i;
+  starts[0]=0;
+  // rows
+  int numberRows = coinModel->numberRows();
+  // convexity
+  rowLower[0]=1.0;
+  rowUpper[0]=1.0;
+  convexity_ = numberRows;
+  starts[1]=0;
+  // x
+  rowLower[1]=0.0;
+  rowUpper[1]=0.0;
+  index[0]=xColumn_;
+  element[0]=-1.0;
+  xRow_ = numberRows+1;
+  starts[2]=1;
+  int nAdd=2;
+  if (xColumn_!=yColumn_) {
+    rowLower[2]=0.0;
+    rowUpper[2]=0.0;
+    index[1]=yColumn;
+    element[1]=-1.0;
+    nAdd=3;
+    yRow_ = numberRows+2;
+    starts[3]=2;
+  } else {
+    yRow_=-1;
+    branchingStrategy_=1;
+  }
+  // may be objective
+  assert (xyRow_>=-1);
+  for (i=0;i<nAdd;i++) {
+    CoinBigIndex iStart = starts[i];
+    coinModel->addRow(starts[i+1]-iStart,index+iStart,element+iStart,rowLower[i],rowUpper[i]);
+  }
+  int n=0;
+  // order is LxLy, LxUy, UxLy and UxUy
+  firstLambda_ = coinModel->numberColumns();
+  // bit sloppy as theoretically could be infeasible but otherwise need to do more work
+  double xB[2];
+  double yB[2];
+  const double * lower = coinModel->columnLowerArray();
+  const double * upper = coinModel->columnUpperArray();
+  xB[0]=lower[xColumn_];
+  xB[1]=upper[xColumn_];
+  yB[0]=lower[yColumn_];
+  yB[1]=upper[yColumn_];
+  if (xMeshSize_!=floor(xMeshSize_)) {
+    // not integral
+    xSatisfied_ = CoinMax(xSatisfied_,0.51*xMeshSize_);
+    if (!yMeshSize_) {
+      xySatisfied_ = CoinMax(xySatisfied_,xSatisfied_*CoinMax(fabs(yB[0]),fabs(yB[1])));
+    }
+  }
+  if (yMeshSize_!=floor(yMeshSize_)) {
+    // not integral
+    ySatisfied_ = CoinMax(ySatisfied_,0.51*yMeshSize_);
+    if (!xMeshSize_) {
+      xySatisfied_ = CoinMax(xySatisfied_,ySatisfied_*CoinMax(fabs(xB[0]),fabs(xB[1])));
+    }
+  }
+  // adjust
+  double distance;
+  double steps;
+  if (xMeshSize_) {
+    distance = xB[1]-xB[0];
+    steps = floor ((distance+0.5*xMeshSize_)/xMeshSize_);
+    distance = xB[0]+xMeshSize_*steps;
+    if (fabs(xB[1]-distance)>xSatisfied_) {
+      printf("bad x mesh %g %g %g -> %g\n",xB[0],xMeshSize_,xB[1],distance);
+      //double newValue = CoinMax(fabs(xB[1]-distance),xMeshSize_);
+      //printf("xSatisfied increased to %g\n",newValue);
+      //xSatisfied_ = newValue;
+      //xB[1]=distance;
+      //coinModel->setColUpper(xColumn_,distance);
+    }
+  }
+  if (yMeshSize_) {
+    distance = yB[1]-yB[0];
+    steps = floor ((distance+0.5*yMeshSize_)/yMeshSize_);
+    distance = yB[0]+yMeshSize_*steps;
+    if (fabs(yB[1]-distance)>ySatisfied_) {
+      printf("bad y mesh %g %g %g -> %g\n",yB[0],yMeshSize_,yB[1],distance);
+      //double newValue = CoinMax(fabs(yB[1]-distance),yMeshSize_);
+      //printf("ySatisfied increased to %g\n",newValue);
+      //ySatisfied_ = newValue;
+      //yB[1]=distance;
+      //coinModel->setColUpper(yColumn_,distance);
+    }
+  }
+  for (i=0;i<4;i++) {
+    double x = (i<2) ? xB[0] : xB[1];
+    double y = ((i&1)==0) ? yB[0] : yB[1];
+    columnLower[i]=0.0;
+    columnUpper[i]=2.0;
+    objective[i]=0.0;
+    double value;
+    // xy
+    value=coefficient_*x*y;
+    if (xyRow_>=0) { 
+      if (fabs(value)<1.0e-19)
+        value = 1.0e-19;
+      element[n]=value;
+      index[n++]=xyRow_;
+    } else {
+      objective[i]=value;
+    }
+    // convexity
+    value=1.0;
+    element[n]=value;
+    index[n++]=0+numberRows;
+    // x
+    value=x;
+    if (fabs(value)<1.0e-19)
+      value = 1.0e-19;
+    element[n]=value;
+    index[n++]=1+numberRows;
+    if (xColumn_!=yColumn_) {
+      // y
+      value=y;
+      if (fabs(value)<1.0e-19)
+      value = 1.0e-19;
+      element[n]=value;
+      index[n++]=2+numberRows;
+    }
+    starts[i+1]=n;
+  }
+  for (i=0;i<4;i++) {
+    CoinBigIndex iStart = starts[i];
+    coinModel->addColumn(starts[i+1]-iStart,index+iStart,element+iStart,columnLower[i],
+                         columnUpper[i],objective[i]);
+  }
+  // At least one has to have a mesh
+  if (!xMeshSize_&&(!yMeshSize_||yRow_<0)) {
+    printf("one of x and y must have a mesh size\n");
+    abort();
+  } else if (yRow_>=0) {
+    if (!xMeshSize_)
+      branchingStrategy_ = 2;
+    else if (!yMeshSize_)
+      branchingStrategy_ = 1;
+  }
+  // Now add constraints to link in x and or y to existing ones.
+  bool xDone=false;
+  bool yDone=false;
+  // order is LxLy, LxUy, UxLy and UxUy
+  for (i=numberExistingObjects-1;i>=0;i--) {
+    const OsiObject * obj = objects[i];
+    const OsiBiLinear * obj2 =
+      dynamic_cast <const OsiBiLinear *>(obj) ;
+    if (obj2) {
+      if (xColumn_==obj2->xColumn_&&!xDone) {
+        // make sure y equal
+        double rhs=0.0;
+        CoinBigIndex starts[2];
+        int index[4];
+        double element[4]= {1.0,1.0,-1.0,-1.0};
+        starts[0]=0;
+        starts[1]=4;
+        index[0]=firstLambda_+0;
+        index[1]=firstLambda_+1;
+        index[2]=obj2->firstLambda_+0;
+        index[3]=obj2->firstLambda_+1;
+        coinModel->addRow(4,index,element,rhs,rhs);
+        xDone=true;
+      }
+      if (yColumn_==obj2->yColumn_&&yRow_>=0&&!yDone) {
+        // make sure x equal
+        double rhs=0.0;
+        CoinBigIndex starts[2];
+        int index[4];
+        double element[4]= {1.0,1.0,-1.0,-1.0};
+        starts[0]=0;
+        starts[1]=4;
+        index[0]=firstLambda_+0;
+        index[1]=firstLambda_+2;
+        index[2]=obj2->firstLambda_+0;
+        index[3]=obj2->firstLambda_+2;
+        coinModel->addRow(4,index,element,rhs,rhs);
+        yDone=true;
+      }
+    }
+  }
+}
 
 // Copy constructor 
@@ -4539 +5099 @@
    xyRow_(rhs.xyRow_),
    convexity_(rhs.convexity_),
+   numberExtraRows_(rhs.numberExtraRows_),
+   multiplier_(NULL),
+   extraRow_(NULL),
    chosen_(rhs.chosen_)
 {
+  if (numberExtraRows_) {
+    multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
+    extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
+  }
 }
 
@@ -4574 +5141 @@
     xyRow_ = rhs.xyRow_;
     convexity_ = rhs.convexity_;
+    numberExtraRows_ = rhs.numberExtraRows_;
+    delete [] multiplier_;
+    delete [] extraRow_;
+    if (numberExtraRows_) {
+      multiplier_ = CoinCopyOfArray(rhs.multiplier_,numberExtraRows_);
+      extraRow_ = CoinCopyOfArray(rhs.extraRow_,numberExtraRows_);
+    } else {
+      multiplier_ = NULL;
+      extraRow_ = NULL;
+    }
     chosen_ = rhs.chosen_;
   }
@@ -4582 +5159 @@
 OsiBiLinear::~OsiBiLinear ()
 {
-}
-
+  delete [] multiplier_;
+  delete [] extraRow_;
+}
+// Adds in data for extra row with variable coefficients
+void 
+OsiBiLinear::addExtraRow(int row, double multiplier)
+{
+  int * tempI = new int [numberExtraRows_+1];
+  double * tempD = new double [numberExtraRows_+1];
+  memcpy(tempI,extraRow_,numberExtraRows_*sizeof(int));
+  memcpy(tempD,multiplier_,numberExtraRows_*sizeof(double));
+  tempI[numberExtraRows_]=row;
+  tempD[numberExtraRows_]=multiplier;
+  if (numberExtraRows_)
+    assert (row>tempI[numberExtraRows_-1]);
+  numberExtraRows_++;
+  delete [] extraRow_;
+  extraRow_ = tempI;
+  delete [] multiplier_;
+  multiplier_ = tempD;
+}
+static bool testCoarse=true;
 // Infeasibility - large is 0.5
 double 
@@ -4657 +5254 @@
   double steps;
   bool xSatisfied;
-  double xNew;
+  double xNew=xB[0];
   if (xMeshSize_) {
     if (x<0.5*(xB[0]+xB[1])) {
@@ -4673 +5270 @@
     }
     // but if first coarse grid then only if gap small
-    if (false&&(branchingStrategy_&8)!=0&&xSatisfied&&
+    if (testCoarse&&(branchingStrategy_&8)!=0&&xSatisfied&&
         xB[1]-xB[0]>=xMeshSize_) {
       // but allow if fine grid would allow
@@ -4687 +5284 @@
   }
   bool ySatisfied;
-  double yNew;
+  double yNew=yB[0];
   if (yMeshSize_) {
     if (y<0.5*(yB[0]+yB[1])) {
@@ -4703 +5300 @@
     }
     // but if first coarse grid then only if gap small
-    if (false&&(branchingStrategy_&8)!=0&&ySatisfied&&
+    if (testCoarse&&(branchingStrategy_&8)!=0&&ySatisfied&&
         yB[1]-yB[0]>=yMeshSize_) {
       // but allow if fine grid would allow
@@ -4765 +5362 @@
     xyLambda /= coefficient_;
   }
+  if (0) {
+    // only true with positive values
+    // see if all convexification constraints OK with true
+    assert (xyTrue+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
+    assert (xyTrue+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
+    assert (xyTrue-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
+    assert (xyTrue-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
+    // see if all convexification constraints OK with lambda version
+#if 1
+    assert (xyLambda+1.0e-5>xB[0]*y+yB[0]*x - xB[0]*yB[0]);
+    assert (xyLambda+1.0e-5>xB[1]*y+yB[1]*x - xB[1]*yB[1]);
+    assert (xyLambda-1.0e-5<xB[1]*y+yB[0]*x - xB[1]*yB[0]);
+    assert (xyLambda-1.0e-5<xB[0]*y+yB[1]*x - xB[0]*yB[1]);
+#endif
+    // see if other bound stuff true
+    assert (xyLambda+1.0e-5>xB[0]*y);
+    assert (xyLambda+1.0e-5>yB[0]*x);
+    assert (xyLambda-1.0e-5<xB[1]*y);
+    assert (xyLambda-1.0e-5<yB[1]*x);
+#define SIZE 2
+    if (yColumn_==xColumn_+SIZE) {
+#if SIZE==6
+      double bMax = 2200.0;
+      double bMin = bMax - 100.0;
+      double b[] = {330.0,360.0,380.0,430.0,490.0,530.0};
+#elif SIZE==2
+      double bMax = 1900.0;
+      double bMin = bMax - 200.0;
+      double b[] = {460.0,570.0};
+#else
+      abort();
+#endif
+      double sum =0.0;
+      double sum2 =0.0;
+      int m=xColumn_;
+      double x = info->solution_[m];
+      double xB[2];
+      double yB[2];
+      xB[0]=info->lower_[m];
+      xB[1]=info->upper_[m];
+      for (int i=0;i<SIZE*SIZE;i+=SIZE) {
+        int n = i+SIZE+m;
+        double y = info->solution_[n];
+        yB[0]=info->lower_[n];
+        yB[1]=info->upper_[n];
+        int firstLambda=SIZE*SIZE+2*SIZE+4*i+4*m;
+        double xyLambda=0.0;
+        for (int j=0;j<4;j++) {
+          int iX = j>>1;
+          int iY = j&1;
+          xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
+        }
+        sum += xyLambda*b[i/SIZE];
+        double xyTrue = x*y;
+        sum2 += xyTrue*b[i/SIZE];
+      }
+      if (sum>bMax*x+1.0e-5||sum<bMin*x-1.0e-5) {
+        //if (sum<bMax*x+1.0e-5&&sum>bMin*x-1.0e-5) {
+        printf("bmin*x %g b*w %g bmax*x %g (true) %g\n",bMin*x,sum,bMax*x,sum2);
+        printf("m %d lb %g value %g up %g\n",
+               m,xB[0],x,xB[1]);
+        sum=0.0;
+        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
+          int n = i+SIZE+m;
+          double y = info->solution_[n];
+          yB[0]=info->lower_[n];
+          yB[1]=info->upper_[n];
+          printf("n %d lb %g value %g up %g\n",
+                 n,yB[0],y,yB[1]);
+          int firstLambda=SIZE*SIZE+2*SIZE+4*i+m*4;
+          double xyLambda=0.0;
+          for (int j=0;j<4;j++) {
+            int iX = j>>1;
+            int iY = j&1;
+            xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda+j];
+            printf("j %d l %d new xylambda %g ",j,firstLambda+j,xyLambda);
+          }
+          sum += xyLambda*b[i/SIZE];
+          printf(" - sum now %g\n",sum);
+        }
+      }
+      if (sum2>bMax*x+1.0e-5||sum2<bMin*x-1.0e-5) {
+        printf("bmin*x %g b*x*y %g bmax*x %g (estimate) %g\n",bMin*x,sum2,bMax*x,sum);
+        printf("m %d lb %g value %g up %g\n",
+               m,xB[0],x,xB[1]);
+        sum2=0.0;
+        for (int i=0;i<SIZE*SIZE;i+=SIZE) {
+          int n = i+SIZE+m;
+          double y = info->solution_[n];
+          yB[0]=info->lower_[n];
+          yB[1]=info->upper_[n];
+          printf("n %d lb %g value %g up %g\n",
+                 n,yB[0],y,yB[1]);
+          double xyTrue = x*y;
+          sum2 += xyTrue*b[i/SIZE];
+          printf("xyTrue %g - sum now %g\n",xyTrue,sum2);
+        }
+      }
+    }
+  }
   // If pseudo shadow prices then see what would happen
   //double pseudoEstimate = 0.0;
@@ -4781 +5478 @@
       // == row so move x and y not xy
     }
+  }
+  if ((branchingStrategy_&16)!=0) {
+    // always treat as satisfied!!
+    xSatisfied=true;
+    ySatisfied=true;
+    xyTrue=xyLambda;
   }
   if ( !xSatisfied) {
@@ -4852 +5555 @@
           }
         }
+        if (testCoarse&&(branchingStrategy_&8)!=0&&xB[1]-xB[0]<1.0001*xSatisfied_&&
+            yB[1]-yB[0]<1.0001*ySatisfied_)
+          feasible=true;
         if (feasible) {
           if (xB[1]-xB[0]>=xSatisfied_&&xMeshSize_) {
@@ -4938 +5644 @@
   }
   whichWay=whichWay_;
+  //if (infeasibility_&&priority_==10)
+  //printf("x %d %g %g %g, y %d %g %g %g\n",xColumn_,xB[0],x,xB[1],yColumn_,yB[0],y,yB[1]);
   return infeasibility_;
 }
-
+// Sets infeasibility and other when pseudo shadow prices
+void
+OsiBiLinear::getPseudoShadow(const OsiBranchingInformation * info)
+{
+  // order is LxLy, LxUy, UxLy and UxUy
+  double xB[2];
+  double yB[2];
+  xB[0]=info->lower_[xColumn_];
+  xB[1]=info->upper_[xColumn_];
+  yB[0]=info->lower_[yColumn_];
+  yB[1]=info->upper_[yColumn_];
+  double x = info->solution_[xColumn_];
+  x = CoinMax(x,xB[0]);
+  x = CoinMin(x,xB[1]);
+  double y = info->solution_[yColumn_];
+  y = CoinMax(y,yB[0]);
+  y = CoinMin(y,yB[1]);
+  int j;
+  double xyTrue = x*y;
+  double xyLambda = 0.0;
+  if ((branchingStrategy_&4)==0) {
+    for (j=0;j<4;j++) {
+      int iX = j>>1;
+      int iY = j&1;
+      xyLambda += xB[iX]*yB[iY]*info->solution_[firstLambda_+j];
+    }
+  } else {
+    if (xyRow_>=0) {
+      const double * element = info->elementByColumn_;
+      const int * row = info->row_;
+      const CoinBigIndex * columnStart = info->columnStart_;
+      const int * columnLength = info->columnLength_;
+      for (j=0;j<4;j++) {
+        int iColumn = firstLambda_+j;
+        int iStart = columnStart[iColumn];
+        int iEnd = iStart + columnLength[iColumn];
+        int k=iStart;
+        double sol = info->solution_[iColumn];
+        for (;k<iEnd;k++) {
+          if (xyRow_==row[k])
+            xyLambda += element[k]*sol;
+        }
+      }
+    } else {
+      // objective
+      const double * objective = info->objective_;
+      for (j=0;j<4;j++) {
+        int iColumn = firstLambda_+j;
+        double sol = info->solution_[iColumn];
+        xyLambda += objective[iColumn]*sol;
+      }
+    }
+    xyLambda /= coefficient_;
+  }
+  assert (info->defaultDual_>=0.0);
+  // If we move to xy then we move by coefficient * (xyTrue-xyLambda) on row xyRow_
+  double movement = xyTrue-xyLambda;
+  infeasibility_=0.0;
+  const double * pi = info->pi_;
+  const double * activity = info->rowActivity_;
+  const double * lower = info->rowLower_;
+  const double * upper = info->rowUpper_;
+  double tolerance = info->primalTolerance_;
+  double direction = info->direction_;
+  if (xyRow_>=0) {
+    assert (!boundType_);
+    if (lower[xyRow_]<-1.0e20) 
+      assert (pi[xyRow_]<=1.0e-3);
+    if (upper[xyRow_]>1.0e20) 
+      assert (pi[xyRow_]>=-1.0e-3);
+    double valueP = pi[xyRow_]*direction;
+    // if move makes infeasible then make at least default
+    double newValue = activity[xyRow_] + movement*coefficient_;
+    if (newValue>upper[xyRow_]+tolerance||newValue<lower[xyRow_]-tolerance) 
+      infeasibility_ += fabs(movement*coefficient_)*CoinMax(fabs(valueP),info->defaultDual_);
+  } else {
+    // objective
+    assert (movement>-1.0e-7);
+    infeasibility_ += movement;
+  }
+  for (int i=0;i<numberExtraRows_;i++) {
+    int iRow = extraRow_[i];
+    if (lower[iRow]<-1.0e20) 
+      assert (pi[iRow]<=1.0e-3);
+    if (upper[iRow]>1.0e20) 
+      assert (pi[iRow]>=-1.0e-3);
+    double valueP = pi[iRow]*direction;
+    // if move makes infeasible then make at least default
+    double newValue = activity[iRow] + movement*multiplier_[i];
+    if (newValue>upper[iRow]+tolerance||newValue<lower[iRow]-tolerance) 
+      infeasibility_ += fabs(movement*multiplier_[i])*CoinMax(fabs(valueP),info->defaultDual_);
+  }
+  if (infeasibility_<1.0e-7)
+    infeasibility_=0.0;
+  otherInfeasibility_ = CoinMax(1.0e-12,infeasibility_*10.0);
+}
 // This looks at solution and sets bounds to contain solution
 double
@@ -5476 +6279 @@
   const int * row = matrix->getIndices();
   const CoinBigIndex * columnStart = matrix->getVectorStarts();
-  //const int * columnLength = matrix->getVectorLengths();
+  const int * columnLength = matrix->getVectorLengths();
   // order is LxLy, LxUy, UxLy and UxUy
   double xB[2];
@@ -5497 +6300 @@
     double y = yB[iY];
     CoinBigIndex k = columnStart[j+firstLambda_];
+    CoinBigIndex last = k+columnLength[j+firstLambda_];
     double value;
     // xy
@@ -5531 +6335 @@
       element[k++]=value;
       numberUpdated++;
+    }
+    // Do extra rows
+    for (int i=0;i<numberExtraRows_;i++) {
+      int iRow = extraRow_[i];
+      for (;k<last;k++) {
+        if (row[k]==iRow)
+          break;
+      }
+      assert (k<last);
+      element[k++] = x*y*multiplier_[i];
     }
   }
@@ -6099 +6913 @@
       if (newUp>upper[iRow]+tolerance||newUp<lower[iRow]-tolerance)
         u = CoinMax(u,info->defaultDual_);
-      upEstimate += u*upMovement;
+      upEstimate += u*upMovement*fabs(el2);
       // if down makes infeasible then make at least default
       double newDown = activity[iRow] - downMovement*el2;
       if (newDown>upper[iRow]+tolerance||newDown<lower[iRow]-tolerance)
         d = CoinMax(d,info->defaultDual_);
-      downEstimate += d*downMovement;
+      downEstimate += d*downMovement*fabs(el2);
     }
     if (downEstimate>=upEstimate) {
@@ -6692 +7506 @@
   return nelCreate;
 }
-#endif
 #include "ClpConstraint.hpp"
 #include "ClpConstraintLinear.hpp"
@@ -6705 +7518 @@
                     int mode)
 {
+#ifdef COIN_HAS_ASL
+  // matrix etc will be changed
+  CoinModel coinModel2 = coinModel;
+  if (coinModel2.moreInfo()) {
+    // for now just ampl objective
+    ClpSimplex * model = new ClpSimplex();
+    model->loadProblem(coinModel2);
+    int numberConstraints;
+    ClpConstraint ** constraints=NULL;
+    int type = model->loadNonLinear(coinModel2.moreInfo(),
+                                    numberConstraints,constraints);
+    if (type==1||type==3) {
+      int returnCode = model->nonlinearSLP(numberPasses,deltaTolerance);
+      assert (!returnCode);
+    } else if (type==2||type==4) {
+      int returnCode = model->nonlinearSLP(numberConstraints,constraints,
+                                           numberPasses,deltaTolerance);
+      assert (!returnCode);
+    } else {
+      printf("error or linear - fix %d\n",type);
+    }
+    //exit(66);
+    return model;
+  }
   // first check and set up arrays
   int numberColumns = coinModel.numberColumns();
   int numberRows = coinModel.numberRows();
   // List of nonlinear rows
-  int * which = new int[numberRows+1];
+  int * which = new int[numberRows];
   bool testLinear=false;
   int numberConstraints=0;
   int iColumn;
-  // matrix etc will be changed
-  CoinModel coinModel2 = coinModel;
   bool linearObjective=true;
   int maximumQuadraticElements=0;
@@ -6749 +7584 @@
     for (iColumn=0;iColumn<numberColumns;iColumn++) 
       coinModel2.setObjective(iColumn,0.0);
-    which[numberConstraints++]=-1;
   }
   int iRow;
@@ -6798 +7632 @@
   ClpSimplex * model = new ClpSimplex();
   // return if nothing
-  if (!numberConstraints) {
+  if (!numberConstraints&&linearObjective) {
     delete [] which;
     model->loadProblem(coinModel);
@@ -6814 +7648 @@
   int saveNumber=numberConstraints;
   numberConstraints=0;
+  ClpQuadraticObjective * quadObj = NULL;
   if (!linearObjective) {
     int numberQuadratic=0;
-    int largestColumn=-1;
     CoinZeroN(linearTerm,numberColumns);
     for (iColumn=0;iColumn<numberColumns;iColumn++) {
@@ -6823 +7657 @@
       const char * expr = coinModel.getColumnObjectiveAsString(iColumn);
       if (strcmp(expr,"Numeric")) {
-        largestColumn = CoinMax(largestColumn,iColumn);
         // value*x*y
         char temp[20000];
@@ -6835 +7668 @@
           if (jColumn>=0) {
             columnQuadratic[numberQuadratic]=jColumn;
-            elementQuadratic[numberQuadratic++]=2.0*value; // convention
-            largestColumn = CoinMax(largestColumn,jColumn);
+            if (jColumn!=iColumn)
+              elementQuadratic[numberQuadratic++]=1.0*value; // convention
+            else if (jColumn==iColumn)
+              elementQuadratic[numberQuadratic++]=2.0*value; // convention
           } else if (jColumn==-2) {
             linearTerm[iColumn] = value;
-            // and put in as row -1
-            columnQuadratic[numberQuadratic]=-1;
-            elementQuadratic[numberQuadratic++]=value;
-            largestColumn = CoinMax(largestColumn,iColumn);
           } else {
             printf("bad nonlinear term %s\n",temp);
@@ -6852 +7683 @@
         // linear part
         linearTerm[iColumn] = coinModel.getColumnObjective(iColumn);
-        // and put in as row -1
-        columnQuadratic[numberQuadratic]=-1;
-        elementQuadratic[numberQuadratic++]=linearTerm[iColumn];
-        if (linearTerm[iColumn])
-          largestColumn = CoinMax(largestColumn,iColumn);
       }
     }
     startQuadratic[numberColumns] = numberQuadratic;
-    // here we create ClpConstraint
-    constraints[numberConstraints++] = new ClpConstraintQuadratic(-1, largestColumn+1, numberColumns,
-                                                                  startQuadratic,columnQuadratic,elementQuadratic);
+    quadObj = new ClpQuadraticObjective(linearTerm, numberColumns,
+                                           startQuadratic,columnQuadratic,elementQuadratic);
   }
   int iConstraint;
@@ -6949 +7774 @@
   delete [] which;
   model->loadProblem(coinModel2);
-  int returnCode = model->nonlinearSLP(numberConstraints, constraints,
-                                       numberPasses,deltaTolerance);
-  // See if any integers
-  for (iConstraint=0;iConstraint<saveNumber;iConstraint++) 
-    delete constraints[iConstraint];
-  
+  if (quadObj)
+    model->setObjective(quadObj);
+  delete quadObj;
+  int returnCode;
+  if (numberConstraints) {
+    returnCode = model->nonlinearSLP(numberConstraints, constraints,
+                                     numberPasses,deltaTolerance);
+    for (iConstraint=0;iConstraint<saveNumber;iConstraint++) 
+      delete constraints[iConstraint];
+  } else {
+    returnCode = model->nonlinearSLP(numberPasses,deltaTolerance);
+  }
   delete [] constraints;
   assert (!returnCode);
   return model;
-}
+#else
+  printf("loadNonLinear needs ampl\n");
+  abort();
+  return NULL;
+#endif
+}
+OsiChooseStrongSubset::OsiChooseStrongSubset() :
+  OsiChooseStrong(),
+  numberObjectsToUse_(0)
+{
+}
+
+OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiSolverInterface * solver) :
+  OsiChooseStrong(solver),
+  numberObjectsToUse_(-1)
+{
+}
+
+OsiChooseStrongSubset::OsiChooseStrongSubset(const OsiChooseStrongSubset & rhs) 
+  : OsiChooseStrong(rhs)
+{  
+  numberObjectsToUse_ = -1;
+}
+
+OsiChooseStrongSubset &
+OsiChooseStrongSubset::operator=(const OsiChooseStrongSubset & rhs)
+{
+  if (this != &rhs) {
+    OsiChooseStrong::operator=(rhs);
+    numberObjectsToUse_ = -1;
+  }
+  return *this;
+}
+
+
+OsiChooseStrongSubset::~OsiChooseStrongSubset ()
+{
+}
+
+// Clone
+OsiChooseVariable *
+OsiChooseStrongSubset::clone() const
+{
+  return new OsiChooseStrongSubset(*this);
+}
+// Initialize
+int 
+OsiChooseStrongSubset::setupList ( OsiBranchingInformation *info, bool initialize)
+{
+  assert (solver_==info->solver_);
+  // Only has to work with Clp
+  OsiSolverInterface * solverA = const_cast<OsiSolverInterface *> (solver_);
+  OsiSolverLink * solver = dynamic_cast<OsiSolverLink *> (solverA);
+  assert (solver);
+  int numberObjects = solver->numberObjects();
+  if (numberObjects>numberObjects_) {
+    // redo useful arrays
+    delete [] upTotalChange_;
+    delete [] downTotalChange_;
+    delete [] upNumber_;
+    delete [] downNumber_;
+    numberObjects_ = solver->numberObjects();
+    upTotalChange_ = new double [numberObjects_];
+    downTotalChange_ = new double [numberObjects_];
+    upNumber_ = new int [numberObjects_];
+    downNumber_ = new int [numberObjects_];
+    CoinZeroN(upTotalChange_,numberObjects_);
+    CoinZeroN(downTotalChange_,numberObjects_);
+    CoinZeroN(upNumber_,numberObjects_);
+    CoinZeroN(downNumber_,numberObjects_);
+  }
+  if (numberObjectsToUse_<0) {
+    // Sort objects so bilinear at end
+    OsiObject ** sorted = new OsiObject * [numberObjects];
+    OsiObject ** objects = solver->objects();
+    numberObjects_=0;
+    int numberBiLinear=0;
+    int i;
+    for (i=0;i<numberObjects;i++) {
+      OsiObject * obj = objects[i];
+      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
+      if (!objB)
+        objects[numberObjects_++]=obj;
+      else
+        sorted[numberBiLinear++]=obj;
+    }
+    numberObjectsToUse_ = numberObjects_;
+    for (i=0;i<numberBiLinear;i++) 
+      objects[numberObjects_++]=sorted[i];
+    delete [] sorted;
+    // See if any master objects
+    for (i=0;i<numberObjectsToUse_;i++) {
+      OsiUsesBiLinear * obj = dynamic_cast<OsiUsesBiLinear *> (objects[i]);
+      if (obj)
+        obj->addBiLinearObjects(solver);
+    }
+  }
+  solver->setNumberObjects(numberObjectsToUse_);
+  numberObjects_=numberObjectsToUse_;
+  // Use shadow prices
+  info->defaultDual_=0.0;
+  int numberUnsatisfied=OsiChooseStrong::setupList ( info, initialize);
+  solver->setNumberObjects(numberObjects);
+  numberObjects_=numberObjects;
+  return numberUnsatisfied;
+}
+/* Choose a variable
+   Returns - 
+   -1 Node is infeasible
+   0  Normal termination - we have a candidate
+   1  All looks satisfied - no candidate
+   2  We can change the bound on a variable - but we also have a strong branching candidate
+   3  We can change the bound on a variable - but we have a non-strong branching candidate
+   4  We can change the bound on a variable - no other candidates
+   We can pick up branch from whichObject() and whichWay()
+   We can pick up a forced branch (can change bound) from whichForcedObject() and whichForcedWay()
+   If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
+*/
+int 
+OsiChooseStrongSubset::chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables)
+{
+  int numberObjects = solver->numberObjects();
+  solver->setNumberObjects(numberObjectsToUse_);
+  numberObjects_=numberObjectsToUse_;
+  // Use shadow prices
+  info->defaultDual_=0.0;
+  int returnCode=OsiChooseStrong::chooseVariable(solver,info,fixVariables);
+  solver->setNumberObjects(numberObjects);
+  numberObjects_=numberObjects;
+  return returnCode;
+}
+/** Default Constructor
+
+  Equivalent to an unspecified binary variable.
+*/
+OsiUsesBiLinear::OsiUsesBiLinear ()
+  : OsiSimpleInteger(),
+    numberBiLinear_(0),
+    type_(0),
+    objects_(NULL)
+{
+}
+
+/** Useful constructor
+
+  Loads actual upper & lower bounds for the specified variable.
+*/
+OsiUsesBiLinear::OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type)
+  : OsiSimpleInteger(solver,iColumn),
+    numberBiLinear_(0),
+    type_(type),
+    objects_(NULL)
+{
+  if (type_) {
+    assert(originalLower_==floor(originalLower_+0.5));
+    assert(originalUpper_==floor(originalUpper_+0.5));
+  }
+}
+
+  
+// Useful constructor - passed solver index and original bounds
+OsiUsesBiLinear::OsiUsesBiLinear ( int iColumn, double lower, double upper, int type)
+  : OsiSimpleInteger(iColumn,lower,upper),
+    numberBiLinear_(0),
+    type_(type),
+    objects_(NULL)
+{
+  if (type_) {
+    assert(originalLower_==floor(originalLower_+0.5));
+    assert(originalUpper_==floor(originalUpper_+0.5));
+  }
+}
+
+// Useful constructor - passed simple integer
+OsiUsesBiLinear::OsiUsesBiLinear ( const OsiSimpleInteger &rhs, int type)
+  : OsiSimpleInteger(rhs),
+    numberBiLinear_(0),
+    type_(type),
+    objects_(NULL)
+{
+  if (type_) {
+    assert(originalLower_==floor(originalLower_+0.5));
+    assert(originalUpper_==floor(originalUpper_+0.5));
+  }
+}
+
+// Copy constructor 
+OsiUsesBiLinear::OsiUsesBiLinear ( const OsiUsesBiLinear & rhs)
+  :OsiSimpleInteger(rhs),
+    numberBiLinear_(0),
+    type_(rhs.type_),
+    objects_(NULL)
+
+{
+}
+
+// Clone
+OsiObject *
+OsiUsesBiLinear::clone() const
+{
+  return new OsiUsesBiLinear(*this);
+}
+
+// Assignment operator 
+OsiUsesBiLinear & 
+OsiUsesBiLinear::operator=( const OsiUsesBiLinear& rhs)
+{
+  if (this!=&rhs) {
+    OsiSimpleInteger::operator=(rhs);
+    delete [] objects_;
+    numberBiLinear_ = 0;
+    type_ = rhs.type_;
+    objects_ = NULL;
+  }
+  return *this;
+}
+
+// Destructor 
+OsiUsesBiLinear::~OsiUsesBiLinear ()
+{
+  delete [] objects_;
+}
+// Infeasibility - large is 0.5
+double 
+OsiUsesBiLinear::infeasibility(const OsiBranchingInformation * info, int & whichWay) const
+{
+  assert (type_==0); // just continuous for now
+  double value = info->solution_[columnNumber_];
+  value = CoinMax(value, info->lower_[columnNumber_]);
+  value = CoinMin(value, info->upper_[columnNumber_]);
+  infeasibility_ = 0.0;
+  for (int i=0;i<numberBiLinear_;i++) {
+    OsiBiLinear * obj = dynamic_cast<OsiBiLinear *> (objects_[i]);
+    assert (obj);
+    obj->getPseudoShadow(info);
+    infeasibility_ += objects_[i]->infeasibility();
+  }
+  bool satisfied=false;
+  whichWay=-1;
+  if (infeasibility_<=info->integerTolerance_) {
+    otherInfeasibility_ = 1.0;
+    satisfied=true;
+    infeasibility_ = 0.0;
+  } else {
+    otherInfeasibility_ = 10.0*infeasibility_;
+    if (value-info->lower_[columnNumber_]>
+        info->upper_[columnNumber_]-value)
+      whichWay=1;
+    else
+      whichWay=-1;
+  }
+  if (preferredWay_>=0&&!satisfied)
+    whichWay = preferredWay_;
+  whichWay_=whichWay;
+  return infeasibility_;
+}
+// Creates a branching object
+OsiBranchingObject * 
+OsiUsesBiLinear::createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const 
+{
+  double value = info->solution_[columnNumber_];
+  value = CoinMax(value, info->lower_[columnNumber_]);
+  value = CoinMin(value, info->upper_[columnNumber_]);
+  assert (info->upper_[columnNumber_]>info->lower_[columnNumber_]);
+  double nearest = floor(value+0.5);
+  double integerTolerance = info->integerTolerance_;
+  if (fabs(value-nearest)<integerTolerance) {
+    // adjust value
+    if (nearest!=info->upper_[columnNumber_])
+      value = nearest+2.0*integerTolerance;
+    else
+      value = nearest-2.0*integerTolerance;
+  }
+  OsiBranchingObject * branch = new OsiIntegerBranchingObject(solver,this,way,
+                                                              value,value,value);
+  return branch;
+}
+// This looks at solution and sets bounds to contain solution
+/** More precisely: it first forces the variable within the existing
+    bounds, and then tightens the bounds to fix the variable at the
+    nearest integer value.
+*/
+double
+OsiUsesBiLinear::feasibleRegion(OsiSolverInterface * solver,
+                                 const OsiBranchingInformation * info) const
+{
+  double value = info->solution_[columnNumber_];
+  double newValue = CoinMax(value, info->lower_[columnNumber_]);
+  newValue = CoinMin(newValue, info->upper_[columnNumber_]);
+  solver->setColLower(columnNumber_,newValue);
+  solver->setColUpper(columnNumber_,newValue);
+  return fabs(value-newValue);
+}
+// Add all bi-linear objects
+void 
+OsiUsesBiLinear::addBiLinearObjects(OsiSolverLink * solver)
+{
+  delete [] objects_;
+  numberBiLinear_=0;
+  OsiObject ** objects = solver->objects();
+  int i;
+  int numberObjects = solver->numberObjects();
+  for (i=0;i<numberObjects;i++) {
+    OsiObject * obj = objects[i];
+    OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
+    if (objB) {
+      if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_)
+        numberBiLinear_++;
+    }
+  }
+  if (numberBiLinear_) {
+    objects_ = new OsiObject * [numberBiLinear_];
+    numberBiLinear_=0;
+    for (i=0;i<numberObjects;i++) {
+      OsiObject * obj = objects[i];
+      OsiBiLinear * objB = dynamic_cast<OsiBiLinear *> (obj);
+      if (objB) {
+        if (objB->xColumn()==columnNumber_||objB->yColumn()==columnNumber_) 
+          objects_[numberBiLinear_++] = obj;;
+      }
+    }
+  } else {
+    objects_=NULL;
+  }
+}

trunk/Cbc/src/CbcLinked.hpp

@@ -8 +8 @@
    CglTemporary
 */
-#ifndef COIN_HAS_LINK
-#ifdef COIN_HAS_ASL
-#define COIN_HAS_LINK
-#endif
-#endif
 #include "CoinModel.hpp"
-#ifdef COIN_HAS_LINK
 #include "OsiClpSolverInterface.hpp"
+#include "OsiChooseVariable.hpp"
+#include "CbcFathom.hpp"
 class CbcModel;
 class CoinPackedMatrix;
@@ -21 +17 @@
 class OsiObject;
 class CglStored;
-//#############################################################################
-
 /**
    
@@ -29 +23 @@
 */
 
-class OsiSolverLink : public OsiClpSolverInterface {
+class OsiSolverLink : public CbcOsiSolver {
   
 public:
@@ -63 +57 @@
       heuristic solution
       Returns solution array
+      mode -
+      0 just get continuous
+      1 round and try normal bab
+      2 use defaultBound_ to bound integer variables near current solution
   */
   double * heuristicSolution(int numberPasses,double deltaTolerance,int mode);
@@ -112 +110 @@
   /// Analyze constraints to see which are convex (quadratic)
   void analyzeObjects();
+  /// Add reformulated bilinear constraints
+  void addTighterConstraints();
   /// Objective value of best solution found internally
   inline double bestObjectiveValue() const
@@ -163 +163 @@
   inline int integerPriority() const
   { return integerPriority_;};
+  /// Objective transfer variable if one
+  inline int objectiveVariable() const
+  { return objectiveVariable_;}
   /// Set biLinear priority
   inline void setBiLinearPriority(int value)
@@ -169 +172 @@
   inline int biLinearPriority() const
   { return biLinearPriority_;};
-  /// Set Cbc Model
-  inline void setCbcModel(CbcModel * model)
-  { cbcModel_=model;};
   /// Return CoinModel
   inline const CoinModel * coinModel() const
@@ -177 +177 @@
   /// Set all biLinear priorities on x-x variables
   void setBiLinearPriorities(int value, double meshSize=1.0);
+  /** Set options and priority on all or some biLinear variables
+      1 - on I-I
+      2 - on I-x
+      4 - on x-x
+      or combinations.
+      -1 means leave (for priority value and strategy value)
+  */
+  void setBranchingStrategyOnVariables(int strategyValue, int priorityValue=-1,
+                                       int mode=7);
   /// Set all mesh sizes on x-x variables
   void setMeshSizes(double value);
@@ -212 +221 @@
   /// Copy of quadratic model if one
   ClpSimplex * quadraticModel_;
-  /// Pointer back to CbcModel
-  CbcModel * cbcModel_;
   /// Number of rows with nonLinearities
   int numberNonLinearRows_;
@@ -236 +243 @@
      1 bit (2) - quadratic only in objective (add OA cuts)
      2 bit (4) - convex
-     4 bit (8) - try adding OA cuts
+     3 bit (8) - try adding OA cuts
+     4 bit (16) - add linearized constraints
   */
   int specialOptions2_;
@@ -681 +689 @@
   */
   OsiBiLinear (OsiSolverInterface * solver, int xColumn,
+               int yColumn, int xyRow, double coefficient,
+               double xMesh, double yMesh,
+               int numberExistingObjects=0,const OsiObject ** objects=NULL );
+  
+  /** Useful constructor - 
+      This Adds in rows and variables to construct valid Linked Ordered Set
+      Adds extra constraints to match other x/y
+      So note not const model
+  */
+  OsiBiLinear (CoinModel * coinModel, int xColumn,
                int yColumn, int xyRow, double coefficient,
                double xMesh, double yMesh,
@@ -793 +811 @@
       next bit
       8 set to say don't use in feasible region
+      next bit
+      16 set to say - Always satisfied !!
   */
   inline int branchingStrategy() const
@@ -820 +840 @@
   /// Compute lambdas (third entry in each .B is current value) (nonzero if bad)
   double computeLambdas(const double xB[3], const double yB[3],const double xybar[4],double lambda[4]) const;
+  /// Adds in data for extra row with variable coefficients
+  void addExtraRow(int row, double multiplier);
+  /// Sets infeasibility and other when pseudo shadow prices
+  void getPseudoShadow(const OsiBranchingInformation * info);
 
 protected:
@@ -857 +881 @@
       next bit
       8 set to say don't use in feasible region
+      next bit
+      16 set to say - Always satisfied !!
   */
   int branchingStrategy_;
@@ -875 +901 @@
   /// Convexity row
   int convexity_;
+  /// Number of extra rows (coefficients to be modified)
+  int numberExtraRows_;
+  /// Multiplier for coefficient on row
+  double * multiplier_;
+  /// Row number
+  int * extraRow_;
   /// Which chosen -1 none, 0 x, 1 y
   mutable short chosen_;
@@ -972 +1004 @@
   int numberPoints_;
 };
-/// Define a single integer class - but one where you kep branching until fixed even if satsified
+/// Define a single integer class - but one where you keep branching until fixed even if satisfied
 
 
@@ -1013 +1045 @@
   /// data
   
+};
+/** Define a single variable class which is involved with OsiBiLinear objects.
+    This is used so can make better decision on where to branch as it can look at 
+    all objects.
+
+    This version sees if it can re-use code from OsiSimpleInteger
+    even if not an integer variable.  If not then need to duplicate code.
+*/
+
+
+class OsiUsesBiLinear : public OsiSimpleInteger {
+
+public:
+
+  /// Default Constructor 
+  OsiUsesBiLinear ();
+
+  /// Useful constructor - passed solver index
+  OsiUsesBiLinear (const OsiSolverInterface * solver, int iColumn, int type);
+  
+  /// Useful constructor - passed solver index and original bounds
+  OsiUsesBiLinear (int iColumn, double lower, double upper, int type);
+  
+  /// Useful constructor - passed simple integer
+  OsiUsesBiLinear (const OsiSimpleInteger & rhs, int type);
+  
+  /// Copy constructor 
+  OsiUsesBiLinear ( const OsiUsesBiLinear & rhs);
+   
+  /// Clone
+  virtual OsiObject * clone() const;
+
+  /// Assignment operator 
+  OsiUsesBiLinear & operator=( const OsiUsesBiLinear& rhs);
+
+  /// Destructor 
+  virtual ~OsiUsesBiLinear ();
+  
+  /// Infeasibility - large is 0.5
+  virtual double infeasibility(const OsiBranchingInformation * info, int & whichWay) const;
+  /** Creates a branching object
+
+    The preferred direction is set by \p way, 0 for down, 1 for up.
+  */
+  virtual OsiBranchingObject * createBranch(OsiSolverInterface * solver, const OsiBranchingInformation * info, int way) const;
+  /** Set bounds to fix the variable at the current value.
+
+    Given an current value, set the lower and upper bounds to fix the
+    variable. Returns amount it had to move variable.
+  */
+  virtual double feasibleRegion(OsiSolverInterface * solver, const OsiBranchingInformation * info) const;
+  /// Add all bi-linear objects
+  void addBiLinearObjects(OsiSolverLink * solver);
+protected:
+  /// data
+  /// Number of bilinear objects (maybe could be more general)
+  int numberBiLinear_;
+  /// Type of variable - 0 continuous, 1 integer
+  int type_;
+  /// Objects
+  OsiObject ** objects_;
+};
+/** This class chooses a variable to branch on
+
+    This is just as OsiChooseStrong but it fakes it so only
+    first so many are looked at in this phase
+
+*/
+
+class OsiChooseStrongSubset  : public OsiChooseStrong {
+ 
+public:
+    
+  /// Default Constructor 
+  OsiChooseStrongSubset ();
+
+  /// Constructor from solver (so we can set up arrays etc)
+  OsiChooseStrongSubset (const OsiSolverInterface * solver);
+
+  /// Copy constructor 
+  OsiChooseStrongSubset (const OsiChooseStrongSubset &);
+   
+  /// Assignment operator 
+  OsiChooseStrongSubset & operator= (const OsiChooseStrongSubset& rhs);
+
+  /// Clone
+  virtual OsiChooseVariable * clone() const;
+
+  /// Destructor 
+  virtual ~OsiChooseStrongSubset ();
+
+  /** Sets up strong list and clears all if initialize is true.
+      Returns number of infeasibilities. 
+      If returns -1 then has worked out node is infeasible!
+  */
+  virtual int setupList ( OsiBranchingInformation *info, bool initialize);
+  /** Choose a variable
+      Returns - 
+     -1 Node is infeasible
+     0  Normal termination - we have a candidate
+     1  All looks satisfied - no candidate
+     2  We can change the bound on a variable - but we also have a strong branching candidate
+     3  We can change the bound on a variable - but we have a non-strong branching candidate
+     4  We can change the bound on a variable - no other candidates
+     We can pick up branch from bestObjectIndex() and bestWhichWay()
+     We can pick up a forced branch (can change bound) from firstForcedObjectIndex() and firstForcedWhichWay()
+     If we have a solution then we can pick up from goodObjectiveValue() and goodSolution()
+     If fixVariables is true then 2,3,4 are all really same as problem changed
+  */
+  virtual int chooseVariable( OsiSolverInterface * solver, OsiBranchingInformation *info, bool fixVariables);
+
+  /// Number of objects to use
+  inline int numberObjectsToUse() const
+  { return numberObjectsToUse_;};
+  /// Set number of objects to use
+  inline void setNumberObjectsToUse(int value)
+  { numberObjectsToUse_ = value;};
+
+protected:
+  // Data
+  /// Number of objects to be used (and set in solver)
+  int numberObjectsToUse_;
 };
 
@@ -1147 +1301 @@
   //@}
 };
-#endif
 class ClpSimplex;
 /** Return an approximate solution to a CoinModel.

trunk/Cbc/src/CbcMessage.cpp

@@ -31 +31 @@
   {CBC_STATUS,10,1,"After %d nodes, %d on tree, %g best solution, best possible %g (%.2f seconds)"},
   {CBC_GAP,11,1,"Exiting as integer gap of %g less than %g or %g%%"},
-  {CBC_ROUNDING,12,1,"Integer solution of %g found by heuristic after %d iterations and %d nodes (%.2f seconds)"},
+  {CBC_ROUNDING,12,1,"Integer solution of %g found by %s after %d iterations and %d nodes (%.2f seconds)"},
   {CBC_TREE_SOL,24,1,"Integer solution of %g found by subtree after %d iterations and %d nodes (%.2f seconds)"},
   {CBC_ROOT,13,1,"At root node, %d cuts changed objective from %g to %g in %d passes"},
@@ -49 +49 @@
   {CBC_START_SUB,28,1,"Starting sub-tree for %s - maximum nodes %d"},
   {CBC_END_SUB,29,1,"Ending sub-tree for %s"},
-  {CBC_HEURISTIC_SOLUTION,30,1,"solution of %g found by %s after %2.f seconds"},
+  {CBC_THREAD_STATS,30,1,"%s%? %d used %d times,  waiting to start %g, %?%g cpu time,%? %g waiting for threads, %? %d locks, %g locked, %g waiting for locks"},
   {CBC_CUTS_STATS,31,1,"%d added rows had average density of %g"},
   {CBC_STRONG_STATS,32,1,"Strong branching done %d times (%d iterations), fathomed %d nodes and fixed %d variables"},
@@ -56 +56 @@
   {CBC_HEURISTICS_OFF,36,1,"Heuristics switched off as %d branching objects are of wrong type"},
   {CBC_STATUS2,37,1,"%d nodes, %d on tree, best %g - possible %g depth %d unsat %d its %d (%.2f seconds)"},
+  {CBC_FPUMP1,38,1,"%s"},
+  {CBC_FPUMP2,39,2,"%s"},
   {CBC_DUMMY_END,999999,0,""}
 };

trunk/Cbc/src/CbcMessage.hpp

@@ -54 +54 @@
   CBC_START_SUB,
   CBC_END_SUB,
-  CBC_HEURISTIC_SOLUTION,
+  CBC_THREAD_STATS,
   CBC_CUTS_STATS,
   CBC_STRONG_STATS,
@@ -61 +61 @@
   CBC_HEURISTICS_OFF,
   CBC_STATUS2,
+  CBC_FPUMP1,
+  CBC_FPUMP2,
   CBC_DUMMY_END
 };

trunk/Cbc/src/CbcModel.cpp

@@ -14 +14 @@
 //#define NODE_LOG
 //#define GLOBAL_CUTS_JUST_POINTERS
+#ifndef CLP_FAST_CODE
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+#endif
 #include <cassert>
 #include <cmath>
@@ -36 +41 @@
 #include "CbcBranchDynamic.hpp"
 #include "CbcHeuristic.hpp"
+#include "CbcHeuristicFPump.hpp"
 #include "CbcModel.hpp"
 #include "CbcTreeLocal.hpp"
@@ -48 +54 @@
 #include "CbcCutGenerator.hpp"
 #include "CbcFeasibilityBase.hpp"
+#include "CbcFathom.hpp"
 // include Probing
 #include "CglProbing.hpp"
@@ -60 +67 @@
 #include "CbcCompareActual.hpp"
 #include "CbcTree.hpp"
+//#define CBC_THREAD
+#ifdef CBC_THREAD
+#include <pthread.h>
+#ifndef CLP_FAST_CODE
+//#define CBC_THREAD_DEBUG 1
+#endif
+#ifdef CBC_THREAD_DEBUG 
+#ifdef NDEBUG
+#undef NDEBUG
+#undef assert
+#         define assert(expression) {                              \
+             if (!(expression)) {                                          \
+                throw CoinError(__STRING(expression), __PRETTY_FUNCTION__, \
+                                "", __FILE__, __LINE__);                   \
+             }                                                             \
+          }
+#endif
+#endif
+// To Pass across to doOneNode
+typedef struct {
+  CbcModel * baseModel;
+  CbcModel * thisModel;
+  CbcNode * node; // filled in every time
+  CbcNode * createdNode; // filled in every time on return
+  pthread_t threadIdOfBase;
+  pthread_mutex_t * mutex; // for locking data
+  pthread_mutex_t * mutex2; // for waking up threads
+  pthread_cond_t * condition2; // for waking up thread
+  int returnCode; // -1 available, 0 busy, 1 finished , 2??
+  double timeLocked;
+  double timeWaitingToLock;
+  double timeWaitingToStart;
+  double timeInThread;
+  int numberTimesLocked;
+  int numberTimesUnlocked;
+  int numberTimesWaitingToStart;
+  int saveStuff[2];
+  struct timespec absTime;
+  bool locked;
+#if CBC_THREAD_DEBUG
+  int threadNumber;
+#endif
+} threadStruct;
+static void * doNodesThread(void * voidInfo);
+static void * doCutsThread(void * voidInfo);
+#endif
 /* Various functions local to CbcModel.cpp */
 
@@ -449 +502 @@
   system held by the solver.
 */
-
 void CbcModel::branchAndBound(int doStatistics) 
 
@@ -461 +513 @@
   strongInfo_[2]=0;
   numberStrongIterations_ = 0;
+#ifndef NDEBUG
+  {
+    int i;
+    int n = solver_->getNumCols();
+    const double *lower = solver_->getColLower() ;
+    const double *upper = solver_->getColUpper() ;
+    for (i=0;i<n;i++) {
+      assert (lower[i]<1.0e10);
+      assert (upper[i]>-1.0e10);
+    }
+    n = solver_->getNumRows();
+    lower = solver_->getRowLower() ;
+    upper = solver_->getRowUpper() ;
+    for (i=0;i<n;i++) {
+      assert (lower[i]<1.0e10);
+      assert (upper[i]>-1.0e10);
+    }
+  }
+#endif
   // original solver (only set if pre-processing)
   OsiSolverInterface * originalSolver=NULL;
@@ -466 +537 @@
   OsiObject ** originalObject = NULL;
   // Set up strategies
+#if 0
+  std::string problemName ;
+  solver_->getStrParam(OsiProbName,problemName) ;
+  if (!strcmp(problemName.c_str(),"PP08A")) solver_->activateRowCutDebugger(problemName.c_str()) ;
+#endif
   if (strategy_) {
     // May do preprocessing
@@ -533 +609 @@
         numberObjects_= numberNewIntegers+numberOldIntegers+numberOldOther+nNonInt;
         object_ = new OsiObject * [numberObjects_];
+        delete [] integerVariable_;
         integerVariable_ = new int [numberNewIntegers+numberOldIntegers];
         /*
@@ -943 +1020 @@
   whichGenerator_ = new int[maximumWhich_] ;
   memset(whichGenerator_,0,maximumWhich_*sizeof(int));
-  int currentNumberCuts = 0 ;
   maximumNumberCuts_ = 0 ;
   currentNumberCuts_ = 0 ;
@@ -1005 +1081 @@
   maximumDepthActual_=0;
   numberDJFixed_=0.0;
+  // Do heuristics
+  {
+    double * newSolution = new double [numberColumns] ;
+    double heuristicValue = getCutoff() ;
+    int found = -1; // no solution found
+
+    currentPassNumber_ = 1; // so root heuristics will run
+    int i;
+    for (i = 0;i<numberHeuristics_;i++) {
+      // see if heuristic will do anything
+      double saveValue = heuristicValue ;
+      int ifSol = heuristic_[i]->solution(heuristicValue,
+                                          newSolution);
+      if (ifSol>0) {
+        // better solution found
+        found = i ;
+        incrementUsed(newSolution);
+        // increment number of solutions so other heuristics can test
+        numberSolutions_++;
+        numberHeuristicSolutions_++;
+      } else {
+        heuristicValue = saveValue ;
+      }
+    }
+    currentPassNumber_ = 0;
+    /*
+      Did any of the heuristics turn up a new solution? Record it before we free
+      the vector.
+    */
+    if (found >= 0) { 
+      // For compiler error on terra cluster!
+      if (found<numberHeuristics_)
+        lastHeuristic_ = heuristic_[found];
+      else
+        lastHeuristic_ = heuristic_[0];
+      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
+      CbcTreeLocal * tree 
+          = dynamic_cast<CbcTreeLocal *> (tree_);
+      if (tree)
+        tree->passInSolution(bestSolution_,heuristicValue);
+    }
+    for (i = 0;i<numberHeuristics_;i++) {
+      // delete FPump
+      CbcHeuristicFPump * pump 
+        = dynamic_cast<CbcHeuristicFPump *> (heuristic_[i]);
+      if (pump) {
+        delete pump;
+        numberHeuristics_ --;
+        for (int j=i;j<numberHeuristics_;j++)
+          heuristic_[j] = heuristic_[j+1];
+      }
+    }
+    delete [] newSolution ;
+  } 
   statistics_ = NULL;
   // Do on switch
@@ -1037 +1167 @@
         feasible = solveWithCuts(cuts, 1,
                                  NULL);
-#if 0
-        currentNumberCuts_ = cuts.sizeRowCuts();
-        if (currentNumberCuts_ >= maximumNumberCuts_) {
-          maximumNumberCuts_ = currentNumberCuts;
-          delete [] addedCuts_;
-          addedCuts_ = new CbcCountRowCut * [maximumNumberCuts_];
-        }
-#endif
       }
     }
@@ -1064 +1186 @@
                            *dblParam_[CbcAllowableFractionGap]);
   if (bestObjective_-bestPossibleObjective_ < testGap && getCutoffIncrement()>=0.0) {
-    if (bestPossibleObjective_<getCutoff()) {
+    if (bestPossibleObjective_<getCutoff())
       stoppedOnGap_ = true ;
-      messageHandler()->message(CBC_GAP,messages())
-        << bestObjective_-bestPossibleObjective_
-        << dblParam_[CbcAllowableGap]
-        << dblParam_[CbcAllowableFractionGap]*100.0
-        << CoinMessageEol ;
-      secondaryStatus_ = 2;
-    }
     feasible = false;
   }
@@ -1127 +1242 @@
       }
     }
-#if 0
-    while (anyAction == -1)
-    {
-      // Set objective value (not so obvious if NLP etc)
-      setObjectiveValue(newNode,NULL);
-      if (numberBeforeTrust_==0 ) {
-        anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ;
-      } else {
-        OsiSolverBranch * branches=NULL;
-        anyAction = newNode->chooseDynamicBranch(this,NULL,branches,numberPassesLeft) ;
-        if (anyAction==-3) 
-          anyAction = newNode->chooseBranch(this,NULL,numberPassesLeft) ; // dynamic did nothing
-      }
-      numberPassesLeft--;
-      if (anyAction == -1)
-      { 
-        //if (solverCharacteristics_->solverType()!=4) 
-          feasible = resolve(NULL,10) != 0 ;
-          //else
-          //feasible = solveWithCuts(cuts,maximumCutPassesAtRoot_,NULL);
-      if (problemFeasibility_->feasible(this,0)<0) {
-        feasible=false; // pretend infeasible
-      }
-      reducedCostFix() ;
-      resolved = true ;
-#       ifdef CBC_DEBUG
-      printf("Resolve (root) as something fixed, Obj value %g %d rows\n",
-             solver_->getObjValue(),
-             solver_->getNumRows()) ;
-#       endif
-      if (!feasible) anyAction = -2 ; }
-      if (anyAction == -2||newNode->objectiveValue() >= cutoff)
-      { delete newNode ;
-        newNode = NULL ;
-        feasible = false ; } }
-#else
-  OsiSolverBranch * branches = NULL;
-  anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
-                           NULL,NULL,NULL,branches);
-  if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
-    if (anyAction != -2) {
-      // zap parent nodeInfo
+    OsiSolverBranch * branches = NULL;
+    anyAction = chooseBranch(newNode, numberPassesLeft, NULL, cuts,resolved,
+                             NULL,NULL,NULL,branches);
+    if (anyAction == -2||newNode->objectiveValue() >= cutoff) {
+      if (anyAction != -2) {
+        // zap parent nodeInfo
 #ifdef COIN_DEVELOP
-      printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
-#endif
-      if (newNode->nodeInfo())
-        newNode->nodeInfo()->nullParent();
-    }
-    delete newNode ;
-    newNode = NULL ;
-    feasible = false ;
-  }
-#endif
+        printf("zapping CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
+#endif
+        if (newNode->nodeInfo())
+          newNode->nodeInfo()->nullParent();
+      }
+      delete newNode ;
+      newNode = NULL ;
+      feasible = false ;
+    }
   }
 /*
@@ -1188 +1266 @@
   assert (!newNode || newNode->objectiveValue() <= cutoff) ;
   // Save address of root node as we don't want to delete it
-  CbcNode * rootNode = newNode;
   // initialize for print out
   int lastDepth=0;
@@ -1204 +1281 @@
   a valid solution for use by setBestSolution().
 */
-  CoinWarmStartBasis *lastws = 0 ;
+  CoinWarmStartBasis *lastws = NULL ;
   if (feasible && newNode->branchingObject())
   { if (resolved)
@@ -1291 +1368 @@
   */
   numberLongStrong_=0;
+  double totalTime = 0.0;
+#ifdef CBC_THREAD
+  CbcNode * createdNode=NULL;
+  CbcModel ** threadModel = NULL;
+  pthread_t * threadId = NULL;
+  int * threadCount = NULL;
+  pthread_mutex_t mutex;
+  pthread_cond_t condition_main;
+  pthread_mutex_t condition_mutex;
+  pthread_mutex_t * mutex2 = NULL;
+  pthread_cond_t * condition2 = NULL;
+  threadStruct * threadInfo = NULL;
+  bool locked=false;
+  int threadStats[6];
+  memset(threadStats,0,sizeof(threadStats));
+  double timeWaiting=0.0;
+  // For now just one model
+  if (numberThreads_) {
+    threadId = new pthread_t [numberThreads_];
+    threadCount = new int [numberThreads_];
+    CoinZeroN(threadCount,numberThreads_);
+    pthread_mutex_init(&mutex,NULL);
+    pthread_cond_init(&condition_main,NULL);
+    pthread_mutex_init(&condition_mutex,NULL);
+    threadModel = new CbcModel * [numberThreads_+1];
+    threadInfo = new threadStruct [numberThreads_+1];
+    mutex2 = new pthread_mutex_t [numberThreads_];
+    condition2 = new pthread_cond_t [numberThreads_];
+    // we don't want a strategy object
+    CbcStrategy * saveStrategy = strategy_;
+    strategy_ = NULL;
+    for (int i=0;i<numberThreads_;i++) {
+      pthread_mutex_init(mutex2+i,NULL);
+      pthread_cond_init(condition2+i,NULL);
+      threadId[i]=0;
+      threadInfo[i].baseModel=this;
+      threadModel[i]=new CbcModel(*this);
+#ifdef COIN_HAS_CLP
+      // Solver may need to know about model
+      CbcModel * thisModel = threadModel[i];
+      CbcOsiSolver * solver =
+              dynamic_cast<CbcOsiSolver *>(thisModel->solver()) ;
+      if (solver)
+        solver->setCbcModel(thisModel);
+#endif
+      mutex_ = (void *) (threadInfo+i);
+      threadModel[i]->moveToModel(this,-1);
+      threadInfo[i].thisModel=threadModel[i];
+      threadInfo[i].node=NULL;
+      threadInfo[i].createdNode=NULL;
+      threadInfo[i].threadIdOfBase=pthread_self();
+      threadInfo[i].mutex=&mutex;
+      threadInfo[i].mutex2=mutex2+i;
+      threadInfo[i].condition2=condition2+i;
+      threadInfo[i].returnCode=-1;
+      threadInfo[i].timeLocked=0.0;
+      threadInfo[i].timeWaitingToLock=0.0;
+      threadInfo[i].timeWaitingToStart=0.0;
+      threadInfo[i].timeInThread=0.0;
+      threadInfo[i].numberTimesLocked=0;
+      threadInfo[i].numberTimesUnlocked=0;
+      threadInfo[i].numberTimesWaitingToStart=0;
+      threadInfo[i].locked=false;
+#if CBC_THREAD_DEBUG
+      threadInfo[i].threadNumber=i+2;
+#endif
+      pthread_create(threadId+i,NULL,doNodesThread,threadInfo+i);
+    }
+    strategy_ = saveStrategy;
+    // Do a partial one for base model
+    threadInfo[numberThreads_].baseModel=this;
+    threadModel[numberThreads_]=this;
+    mutex_ = (void *) (threadInfo+numberThreads_);
+    threadInfo[numberThreads_].node=NULL;
+    threadInfo[numberThreads_].mutex=&mutex;
+    threadInfo[numberThreads_].condition2=&condition_main;
+    threadInfo[numberThreads_].mutex2=&condition_mutex;
+    threadInfo[numberThreads_].timeLocked=0.0;
+    threadInfo[numberThreads_].timeWaitingToLock=0.0;
+    threadInfo[numberThreads_].numberTimesLocked=0;
+    threadInfo[numberThreads_].numberTimesUnlocked=0;
+    threadInfo[numberThreads_].locked=false;
+#if CBC_THREAD_DEBUG
+    threadInfo[numberThreads_].threadNumber=1;
+#endif
+  }
+#endif
 /*
   At last, the actual branch-and-cut search loop, which will iterate until
@@ -1301 +1465 @@
   than the current objective cutoff.
 */
-  while (!tree_->empty()) {
+  while (true) {
+#ifdef CBC_THREAD
+    if (!locked) {
+      lockThread();
+      locked=true;
+    }
+#endif
+    if (tree_->empty()) {
+#ifdef CBC_THREAD
+      if (numberThreads_) {
+#ifdef COIN_DEVELOP
+        printf("empty\n");
+#endif
+        // may still be outstanding nodes
+        int iThread;
+        for (iThread=0;iThread<numberThreads_;iThread++) {
+          if (threadId[iThread]) {
+            if (threadInfo[iThread].returnCode==0) 
+              break;
+          }
+        }
+        if (iThread<numberThreads_) {
+#ifdef COIN_DEVELOP
+          printf("waiting for thread %d code 0\n",iThread);
+#endif
+          unlockThread();
+          locked = false;
+          pthread_cond_signal(threadInfo[iThread].condition2); // unlock in case
+          while (true) {
+            pthread_mutex_lock(&condition_mutex);
+            struct timespec absTime;
+            clock_gettime(CLOCK_REALTIME,&absTime);
+            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
+            absTime.tv_nsec += 1000000; // millisecond
+            if (absTime.tv_nsec>=1000000000) {
+              absTime.tv_nsec -= 1000000000;
+              absTime.tv_sec++;
+            }
+            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
+            clock_gettime(CLOCK_REALTIME,&absTime);
+            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
+            timeWaiting += time2-time;
+            pthread_mutex_unlock(&condition_mutex);
+            if (threadInfo[iThread].returnCode!=0) 
+              break;
+            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+          }
+          threadModel[iThread]->moveToModel(this,1);
+          assert (threadInfo[iThread].returnCode==1);
+          // say available
+          threadInfo[iThread].returnCode=-1;
+          threadStats[4]++;
+#ifdef COIN_DEVELOP
+          printf("thread %d code now -1\n",iThread);
+#endif
+          continue;
+        } else {
+#ifdef COIN_DEVELOP
+          printf("no threads at code 0 \n");
+#endif
+          // now check if any have just finished
+          for (iThread=0;iThread<numberThreads_;iThread++) {
+            if (threadId[iThread]) {
+              if (threadInfo[iThread].returnCode==1) 
+                break;
+            }
+          }
+          if (iThread<numberThreads_) {
+            unlockThread();
+            locked = false;
+            threadModel[iThread]->moveToModel(this,1);
+            assert (threadInfo[iThread].returnCode==1);
+            // say available
+            threadInfo[iThread].returnCode=-1;
+            threadStats[4]++;
+#ifdef COIN_DEVELOP
+            printf("thread %d code now -1\n",iThread);
+#endif
+            continue;
+          }
+        }
+        if (!tree_->empty()) {
+#ifdef COIN_DEVELOP
+          printf("tree not empty!!!!!!\n");
+#endif
+          continue;
+        }
+        for (iThread=0;iThread<numberThreads_;iThread++) {
+          if (threadId[iThread]) {
+            if (threadInfo[iThread].returnCode!=-1) { 
+              printf("bad end of tree\n");
+              abort();
+            }
+          }
+        }
+#ifdef COIN_DEVELOP
+        printf("finished ************\n");
+#endif
+      }
+      unlockThread();
+      locked=false; // not needed as break
+#endif
+      break;
+    }
+#ifdef CBC_THREAD
+    unlockThread();
+    locked = false;
+#endif
+/*
+  Check for abort on limits: node count, solution count, time, integrality gap.
+*/
+    totalTime = getCurrentSeconds() ;
+    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
+          numberSolutions_ < intParam_[CbcMaxNumSol] &&
+          totalTime < dblParam_[CbcMaximumSeconds] &&
+          !stoppedOnGap_&&!eventHappened_)) {
+      // out of loop
+      break;
+    }
 #ifdef BONMIN
     assert(!solverCharacteristics_->solutionAddsCuts() || solverCharacteristics_->mipFeasible());
@@ -1308 +1590 @@
       double newCutoff = getCutoff();
       if (analyzeResults_) {
-        // see if we could fix any (more)
-        int n=0;
-        double * newLower = analyzeResults_;
-        double * objLower = newLower+numberIntegers_;
-        double * newUpper = objLower+numberIntegers_;
-        double * objUpper = newUpper+numberIntegers_;
-        for (int i=0;i<numberIntegers_;i++) {
-          if (objLower[i]>newCutoff) {
-            n++;
-            if (objUpper[i]>newCutoff) {
-              newCutoff = -COIN_DBL_MAX;
-              break;
-            }
-          } else if (objUpper[i]>newCutoff) {
-            n++;
-          }
-        }
-        if (newCutoff==-COIN_DBL_MAX) {
-          printf("Root analysis says finished\n");
-        } else if (n>numberFixedNow_) {
-          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
-          numberFixedNow_=n;
-        }
+        // see if we could fix any (more)
+        int n=0;
+        double * newLower = analyzeResults_;
+        double * objLower = newLower+numberIntegers_;
+        double * newUpper = objLower+numberIntegers_;
+        double * objUpper = newUpper+numberIntegers_;
+        for (int i=0;i<numberIntegers_;i++) {
+          if (objLower[i]>newCutoff) {
+            n++;
+            if (objUpper[i]>newCutoff) {
+              newCutoff = -COIN_DBL_MAX;
+              break;
+            }
+          } else if (objUpper[i]>newCutoff) {
+            n++;
+          }
+        }
+        if (newCutoff==-COIN_DBL_MAX) {
+          printf("Root analysis says finished\n");
+        } else if (n>numberFixedNow_) {
+          printf("%d more fixed by analysis - now %d\n",n-numberFixedNow_,n);
+          numberFixedNow_=n;
+        }
       }
       if (eventHandler) {
-        if (!eventHandler->event(CbcEventHandler::solution)) {
-          eventHappened_=true; // exit
-        }
-      }
+        if (!eventHandler->event(CbcEventHandler::solution)) {
+          eventHappened_=true; // exit
+        }
+      }
+      lockThread();
       // Do from deepest
       tree_->cleanTree(this, newCutoff,bestPossibleObjective_) ;
       nodeCompare_->newSolution(this) ;
       nodeCompare_->newSolution(this,continuousObjective_,
-                                continuousInfeasibilities_) ;
+                                continuousInfeasibilities_) ;
       tree_->setComparison(*nodeCompare_) ;
-      if (tree_->empty())
-        break; // finished
+      if (tree_->empty()) {
+        unlockThread();
+        // For threads we need to check further
+        //break; // finished
+        continue;
+      }
+      unlockThread();
     }
     cutoff = getCutoff() ;
 /*
-  Periodic activities: Opportunities to
+    Periodic activities: Opportunities to
     + tweak the nodeCompare criteria,
     + check if we've closed the integrality gap enough to quit, 
@@ -1354 +1642 @@
 */
     if ((numberNodes_%1000) == 0) {
+      lockThread();
       bool redoTree=nodeCompare_->every1000Nodes(this, numberNodes_) ;
 #ifdef CHECK_CUT_SIZE
@@ -1361 +1650 @@
       if (redoTree)
         tree_->setComparison(*nodeCompare_) ;
+      unlockThread();
     }
     if (saveCompare&&!hotstartSolution_) {
@@ -1368 +1658 @@
       saveCompare=NULL;
       // redo tree
+      lockThread();
       tree_->setComparison(*nodeCompare_) ;
+      unlockThread();
     }
     if ((numberNodes_%printFrequency_) == 0) {
+      lockThread();
       int j ;
       int nNodes = tree_->size() ;
@@ -1379 +1672 @@
           bestPossibleObjective_ = node->objectiveValue() ;
       }
+      unlockThread();
       if (!intParam_[CbcPrinting]) {
         messageHandler()->message(CBC_STATUS,messages())
@@ -1405 +1699 @@
       stoppedOnGap_ = true ;
     }
-
+    
 #   ifdef CHECK_NODE_FULL
     verifyTreeNodes(tree_,*this) ;
@@ -1412 +1706 @@
     verifyCutCounts(tree_,*this) ;
 #   endif
-
 /*
   Now we come to the meat of the loop. To create the active subproblem, we'll
@@ -1423 +1716 @@
     if (!node)
       continue;
+#ifndef CBC_THREAD
+    int currentNumberCuts = 0 ;
     currentNode_=node; // so can be accessed elsewhere
 #ifdef CBC_DEBUG
@@ -1429 +1724 @@
            node->guessedObjectiveValue());
 #endif
+#if NEW_UPDATE_OBJECT==0
     // Save clone in branching decision
     if(branchingMethod_)
       branchingMethod_->saveBranchingObject(node->modifiableBranchingObject());
-    bool nodeOnTree=false; // Node has been popped
+#endif
     // Say not on optimal path
     bool onOptimalPath=false;
@@ -1460 +1756 @@
     CbcNodeInfo * nodeInfo = node->nodeInfo() ;
     newNode = NULL ;
+    int branchesLeft=0;
     if (!addCuts(node,lastws,numberFixedNow_>numberFixedAtRoot_))
     { int i ;
@@ -1471 +1768 @@
         solverCharacteristics_->setBeforeUpper(upperBefore);
       }
-      bool deleteNode ;
       if (messageHandler()->logLevel()>2)
         node->modifiableBranchingObject()->print();
-      int retCode;
       if (!useOsiBranching) 
-        retCode = node->branch(NULL); // old way
+        branchesLeft = node->branch(NULL); // old way
       else
-        retCode = node->branch(solver_); // new way
-      if (retCode)
-      { 
+        branchesLeft = node->branch(solver_); // new way
+      if (branchesLeft) {
         // set nodenumber correctly
         node->nodeInfo()->setNodeNumber(numberNodes2_);
@@ -1497 +1791 @@
         }
         numberNodes2_++;
-        nodeOnTree=true; // back on tree
-        deleteNode = false ;
+        //nodeOnTree=true; // back on tree
+        //deleteNode = false ;
 #       ifdef CHECK_NODE
         printf("Node %x pushed back on tree - %d left, %d count\n",node,
@@ -1504 +1798 @@
                nodeInfo->numberPointingToThis()) ;
 #       endif
-      }
-      else
-      { deleteNode = true ;
-      if (!nodeInfo->numberBranchesLeft())
-        nodeInfo->allBranchesGone(); // can clean up
-      }
-
+      } else {
+        //deleteNode = true ;
+        if (!nodeInfo->numberBranchesLeft())
+          nodeInfo->allBranchesGone(); // can clean up
+      }
       if ((specialOptions_&1)!=0) {
         /*
@@ -1524 +1816 @@
         }
       }
+      
 /*
   Reoptimize, possibly generating cuts and/or using heuristics to find
@@ -1597 +1890 @@
       }
 /*
-  Check for abort on limits: node count, solution count, time, integrality gap.
-*/
-      double totalTime = CoinCpuTime()-dblParam_[CbcStartSeconds] ;
-      if (numberNodes_ < intParam_[CbcMaxNumNode] &&
-          numberSolutions_ < intParam_[CbcMaxNumSol] &&
-          totalTime < dblParam_[CbcMaximumSeconds] &&
-          !stoppedOnGap_&&!eventHappened_) 
-      {
-/*
   Are we still feasible? If so, create a node and do the work to attach a
   branching object, reoptimising as needed if chooseBranch() identifies
@@ -1645 +1929 @@
           int numberPassesLeft=5;
           checkingNode=true;
-#if 0
-          while (anyAction == -1)
-          { 
-            // Set objective value (not so obvious if NLP etc)
-            setObjectiveValue(newNode,node);
-            if (numberBeforeTrust_==0 ) {
-              anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ;
-            } else {
-              OsiSolverBranch * branches=NULL;
-              anyAction = newNode->chooseDynamicBranch(this,node,branches,numberPassesLeft) ;
-              if (anyAction==-3) 
-                anyAction = newNode->chooseBranch(this,node,numberPassesLeft) ; // dynamic did nothing
-            }
-            if (solverCharacteristics_ && 
-                solverCharacteristics_->solutionAddsCuts() && // we are in some OA based bab
-                feasible && (newNode->numberUnsatisfied()==0) //solution has become integer feasible during strong branching
-                )
-              { 
-                //in the present case we need to check here integer infeasibility if the node is not fathomed we will have to do the loop
-                // again
-                std::cout<<solver_<<std::endl;
-                resolve(solver_);
-                double objval = solver_->getObjValue();
-                setBestSolution(CBC_SOLUTION, objval,
-                                solver_->getColSolution()) ;
-                lastHeuristic_ = NULL;
-                int easy=2;
-                if (!solverCharacteristics_->mipFeasible())//did we prove that the node could be pruned?
-                  feasible = false;
-                // Reset the bound now
-                solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
-                
-                
-                //numberPassesLeft++;
-                solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
-                feasible &= resolve(node ? node->nodeInfo() : NULL,11) != 0 ;
-                solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
-                resolved = true ;
-                if (problemFeasibility_->feasible(this,0)<0) {
-                  feasible=false; // pretend infeasible
-                }
-                if(feasible)
-                  anyAction = -1;
-                else
-                  anyAction = -2;
-              }
-/*
-  Yep, false positives for sure. And no easy way to distinguish honest
-  infeasibility from `found a solution and tightened objective target.'
-
-            if (onOptimalPath)
-              assert (anyAction!=-2); // can be useful but gives false positives on strong
-*/
-            numberPassesLeft--;
-            if (numberPassesLeft<=-1) {
-              if (!numberLongStrong_)
-                messageHandler()->message(CBC_WARNING_STRONG,
-                                          messages()) << CoinMessageEol ;
-              numberLongStrong_++;
-            }
-            if (anyAction == -1)
-            {
-              // can do quick optimality check
-              int easy=2;
-              solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,&easy) ;
-              //if (solverCharacteristics_->solverType()!=4) 
-                feasible = resolve(node ? node->nodeInfo() : NULL,11) != 0 ;
-                //else
-                //feasible = solveWithCuts(cuts,maximumCutPasses_,node);
-              solver_->setHintParam(OsiDoInBranchAndCut,true,OsiHintDo,NULL) ;
-              resolved = true ;
-              if (problemFeasibility_->feasible(this,0)<0) {
-                feasible=false; // pretend infeasible
-              }
-              if (feasible)
-              { 
-                // Set objective value (not so obvious if NLP etc)
-                setObjectiveValue(newNode,node);
-                reducedCostFix() ;
-                if (newNode->objectiveValue() >= getCutoff()) 
-                  anyAction=-2;
-              }
-              else
-              { anyAction = -2 ; } } }
-          if (anyAction >= 0)
-          { if (resolved)
-            { bool needValidSolution = (newNode->branchingObject() == NULL) ;
-              takeOffCuts(cuts,needValidSolution,NULL) ; 
-#             ifdef CHECK_CUT_COUNTS
-              { printf("Number of rows after chooseBranch fix (node)"
-                       "(active only) %d\n",
-                        numberRowsAtContinuous_+numberNewCuts_+
-                        numberOldActiveCuts_) ;
-                const CoinWarmStartBasis* debugws =
-                  dynamic_cast<const CoinWarmStartBasis*>
-                    (solver_->getWarmStart()) ;
-                debugws->print() ;
-                delete debugws ; }
-#             endif
-            }
-            newNode->createInfo(this,node,lastws,lowerBefore,upperBefore,
-                                numberOldActiveCuts_,numberNewCuts_) ;
-            if (newNode->numberUnsatisfied()) {
-              maximumDepthActual_ = CoinMax(maximumDepthActual_,newNode->depth());
-              newNode->initializeInfo() ;
-              newNode->nodeInfo()->addCuts(cuts,newNode->numberBranches(),
-                                           whichGenerator_) ; } } }
-        else {
-          anyAction = -2 ; 
-          // Reset bound anyway (no harm if not odd)
-          solverCharacteristics_->setMipBound(-COIN_DBL_MAX);
-          node->nodeInfo()->decrement();
-        }
-        // May have slipped through i.e. anyAction == 0 and objective above cutoff
-        // I think this will screw up cut reference counts if executed.
-        // We executed addCuts just above. (lh)
-        if ( anyAction >=0 ) {
-          assert (newNode);
-          if (newNode->objectiveValue() >= getCutoff()) {
-            anyAction = -2; // say bad after all
-#ifdef COIN_DEVELOP
-            printf("zapping2 CbcNodeInfo %x\n",newNode->nodeInfo()->parent());
-#endif
-            // zap parent nodeInfo
-            if (newNode->nodeInfo())
-              newNode->nodeInfo()->nullParent();
-          }
-        }
-/*
-  If we end up infeasible, we can delete the new node immediately. Since this
-  node won't be needing the cuts we collected, decrement the reference counts.
-  If we are feasible, then we'll be placing this node into the live set, so
-  increment the reference count in the current (parent) nodeInfo.
-*/
-        if (anyAction == -2)
-        { delete newNode ;
-          newNode = NULL ;
-          // say strong doing well
-          if (checkingNode)
-            setSpecialOptions(specialOptions_|8);
-          for (i = 0 ; i < currentNumberCuts_ ; i++)
-          { if (addedCuts_[i])
-            { if (!addedCuts_[i]->decrement(1))
-                delete addedCuts_[i] ; } } }
-        else
-        { nodeInfo->increment() ;
-        if ((numberNodes_%20)==0) {
-          // say strong not doing as well
-          setSpecialOptions(specialOptions_&~8);
-        }
-#else
         OsiSolverBranch * branches=NULL;
         // point to useful information
@@ -1823 +1956 @@
           }
         }
-#endif
         }
 /*
@@ -1866 +1998 @@
                 generate=false; // only special cuts
               if (generate) {
-                generator_[i]->generateCuts(theseCuts,true,NULL) ;
+                generator_[i]->generateCuts(theseCuts,true,solver_,NULL) ;
                 int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
                 if (numberRowCutsAfter) {
@@ -1913 +2045 @@
                 heurValue = saveValue ; } }
             if (found >= 0) {
+              lastHeuristic_ = heuristic_[found];
               setBestSolution(CBC_ROUNDING,heurValue,newSolution) ;
-              lastHeuristic_ = heuristic_[found];
             }
             delete [] newSolution ;
@@ -1964 +2096 @@
   set.
 */
-        if (deleteNode) 
-          delete node ; }
-/*
-  End of the non-abort actions. The next block of code is executed if we've
-  aborted because we hit one of the limits. Clean up by deleting the live set
-  and break out of the node processing loop. Note that on an abort, node may
-  have been pushed back onto the tree for further processing, in which case
-  it'll be deleted in cleanTree. We need to check.
-*/
+      if (branchesLeft)
+      { 
+      }
       else
       { 
-        if (tree_->size())
-          tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
-        delete nextRowCut_;
-        // We need to get rid of node if is has already been popped from tree
-        if (!nodeOnTree&&!stoppedOnGap_&&node!=rootNode)
-          delete node;
-        if (stoppedOnGap_)
-        { messageHandler()->message(CBC_GAP,messages())
-            << bestObjective_-bestPossibleObjective_
-            << dblParam_[CbcAllowableGap]
-            << dblParam_[CbcAllowableFractionGap]*100.0
-            << CoinMessageEol ;
-        secondaryStatus_ = 2;
-          status_ = 0 ; }
-        else
-        if (isNodeLimitReached())
-        { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
-        secondaryStatus_ = 3;
-          status_ = 1 ; }
-        else
-        if (totalTime >= dblParam_[CbcMaximumSeconds])
-        { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
-        secondaryStatus_ = 4;
-          status_ = 1 ; }
-        else
-        if (eventHappened_)
-        { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
-        secondaryStatus_ = 5;
-          status_ = 5 ; }
-        else
-        { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
-        secondaryStatus_ = 6;
-          status_ = 1 ; }
-        break ; }
+        if (!nodeInfo->numberBranchesLeft())
+          nodeInfo->allBranchesGone(); // can clean up
+        delete node ; }
+    } else {
+      // add cuts found to be infeasible (on bound)!
+      abort();
+      delete node;
+    }
 /*
   Delete cuts to get back to the original system.
@@ -2022 +2122 @@
         { delRows[i] = i+numberRowsAtContinuous_ ; }
         solver_->deleteRows(numberToDelete,delRows) ;
-        delete [] delRows ; } }
-/*
-  This node fathomed when addCuts atttempted to revive it. Toss it.
-*/
-    else
-      { delete node ; } }
+        delete [] delRows ; }
+#else // end of not CBC_THREAD
+      if (!numberThreads_) {
+        doOneNode(this,node,createdNode);
+      } else {
+        threadStats[0]++;
+        //need to think
+        int iThread;
+        // Start one off if any available
+        for (iThread=0;iThread<numberThreads_;iThread++) {
+          if (threadInfo[iThread].returnCode==-1) 
+            break;
+        }
+        if (iThread<numberThreads_) {
+          threadInfo[iThread].node=node;
+          assert (threadInfo[iThread].returnCode==-1);
+          // say in use
+          threadInfo[iThread].returnCode=0;
+          threadModel[iThread]->moveToModel(this,0);
+          pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+          threadCount[iThread]++;
+        }
+        lockThread();
+        locked=true;
+        // see if any finished
+        for (iThread=0;iThread<numberThreads_;iThread++) {
+          if (threadInfo[iThread].returnCode>0) 
+            break;
+        }
+        unlockThread();
+        locked=false;
+        if (iThread<numberThreads_) {
+          threadModel[iThread]->moveToModel(this,1);
+          assert (threadInfo[iThread].returnCode==1);
+          // say available
+          threadInfo[iThread].returnCode=-1;
+          // carry on
+          threadStats[3]++;
+        } else {
+          // Start one off if any available
+          for (iThread=0;iThread<numberThreads_;iThread++) {
+            if (threadInfo[iThread].returnCode==-1) 
+              break;
+          }
+          if (iThread<numberThreads_) {
+            lockThread();
+            locked=true;
+            // If any on tree get
+            if (!tree_->empty()) {
+              //node = tree_->bestNode(cutoff) ;
+              //assert (node);
+              threadStats[1]++;
+              continue; // ** get another node
+            }
+            unlockThread();
+            locked=false;
+          }
+          // wait (for debug could sleep and use test)
+          bool finished=false;
+          while (!finished) {
+            pthread_mutex_lock(&condition_mutex);
+            struct timespec absTime;
+            clock_gettime(CLOCK_REALTIME,&absTime);
+            double time = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
+            absTime.tv_nsec += 1000000; // millisecond
+            if (absTime.tv_nsec>=1000000000) {
+              absTime.tv_nsec -= 1000000000;
+              absTime.tv_sec++;
+            }
+            pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
+            clock_gettime(CLOCK_REALTIME,&absTime);
+            double time2 = absTime.tv_sec+1.0e-9*absTime.tv_nsec;
+            timeWaiting += time2-time;
+            pthread_mutex_unlock(&condition_mutex);
+            for (iThread=0;iThread<numberThreads_;iThread++) {
+              if (threadInfo[iThread].returnCode>0) {
+                finished=true;
+                break;
+              } else if (threadInfo[iThread].returnCode==0) {
+                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+              }
+            }
+          }
+          assert (iThread<numberThreads_);
+          threadModel[iThread]->moveToModel(this,1);
+          node = threadInfo[iThread].node;
+          threadInfo[iThread].node=NULL;
+          assert (threadInfo[iThread].returnCode==1);
+          // say available
+          threadInfo[iThread].returnCode=-1;
+          // carry on
+          threadStats[2]++;
+        }
+      }
+      //lastDepth=node->depth();
+      //lastUnsatisfied=node->numberUnsatisfied();
+#endif // end of CBC_THREAD
+  }
+#ifdef CBC_THREAD
+  if (numberThreads_) {
+    //printf("stats ");
+    //for (unsigned int j=0;j<sizeof(threadStats)/sizeof(int);j++)
+    //printf("%d ",threadStats[j]);
+    //printf("\n");
+    int i;
+    // Seems to be bug in CoinCpu on Linux - does threads as well despite documentation
+    double time=0.0;
+    for (i=0;i<numberThreads_;i++) 
+      time += threadInfo[i].timeInThread;
+    bool goodTimer = time<(getCurrentSeconds());
+    //bool stopped = (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
+    //        numberSolutions_ < intParam_[CbcMaxNumSol] &&
+    //        totalTime < dblParam_[CbcMaximumSeconds] &&
+    //        !stoppedOnGap_&&!eventHappened_));
+    for (i=0;i<numberThreads_;i++) {
+      while (threadInfo[i].returnCode==0) {
+        pthread_cond_signal(threadInfo[i].condition2); // unlock
+        pthread_mutex_lock(&condition_mutex);
+        struct timespec absTime;
+        clock_gettime(CLOCK_REALTIME,&absTime);
+        absTime.tv_nsec += 1000000; // millisecond
+        if (absTime.tv_nsec>=1000000000) {
+          absTime.tv_nsec -= 1000000000;
+          absTime.tv_sec++;
+        }
+        pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
+        clock_gettime(CLOCK_REALTIME,&absTime);
+        pthread_mutex_unlock(&condition_mutex);
+      }
+      pthread_cond_signal(threadInfo[i].condition2); // unlock
+      pthread_mutex_lock(&condition_mutex); // not sure necessary but have had one hang on interrupt
+      threadModel[i]->numberThreads_=0; // say exit
+      threadInfo[i].returnCode=0;
+      pthread_mutex_unlock(&condition_mutex);
+      pthread_cond_signal(threadInfo[i].condition2); // unlock
+      //if (!stopped)
+        pthread_join(threadId[i],NULL);
+        //else
+        //pthread_kill(threadId[i]); // kill rather than try and synchronize
+      threadModel[i]->moveToModel(this,2);
+      pthread_mutex_destroy (threadInfo[i].mutex2);
+      pthread_cond_destroy (threadInfo[i].condition2);
+      assert (threadInfo[i].numberTimesLocked==threadInfo[i].numberTimesUnlocked);
+      handler_->message(CBC_THREAD_STATS,messages_)
+        <<"Thread";
+      handler_->printing(true)
+        <<i<<threadCount[i]<<threadInfo[i].timeWaitingToStart;
+      handler_->printing(goodTimer)<<threadInfo[i].timeInThread;
+      handler_->printing(false)<<0.0;
+      handler_->printing(true)<<threadInfo[i].numberTimesLocked
+        <<threadInfo[i].timeLocked<<threadInfo[i].timeWaitingToLock
+        <<CoinMessageEol;
+    }
+    assert (threadInfo[numberThreads_].numberTimesLocked==threadInfo[numberThreads_].numberTimesUnlocked);
+    handler_->message(CBC_THREAD_STATS,messages_)
+      <<"Main thread";
+    handler_->printing(false)<<0<<0<<0.0;
+    handler_->printing(false)<<0.0;
+    handler_->printing(true)<<timeWaiting;
+    handler_->printing(true)<<threadInfo[numberThreads_].numberTimesLocked
+      <<threadInfo[numberThreads_].timeLocked<<threadInfo[numberThreads_].timeWaitingToLock
+      <<CoinMessageEol;
+    pthread_mutex_destroy (&mutex);
+    pthread_cond_destroy (&condition_main);
+    pthread_mutex_destroy (&condition_mutex);
+    // delete models (here in case some point to others)
+    for (i=0;i<numberThreads_;i++) {
+      delete threadModel[i];
+    }
+    delete [] mutex2;
+    delete [] condition2;
+    delete [] threadId;
+    delete [] threadInfo;
+    delete [] threadModel;
+    delete [] threadCount;
+    mutex_=NULL;
+    // adjust time to allow for children on some systems
+    dblParam_[CbcStartSeconds] -= CoinCpuTimeJustChildren();
+  }
+#endif
+/*
+  End of the non-abort actions. The next block of code is executed if we've
+  aborted because we hit one of the limits. Clean up by deleting the live set
+  and break out of the node processing loop. Note that on an abort, node may
+  have been pushed back onto the tree for further processing, in which case
+  it'll be deleted in cleanTree. We need to check.
+*/
+    if (!(numberNodes_ < intParam_[CbcMaxNumNode] &&
+        numberSolutions_ < intParam_[CbcMaxNumSol] &&
+        totalTime < dblParam_[CbcMaximumSeconds] &&
+        !stoppedOnGap_&&!eventHappened_)) {
+      if (tree_->size())
+        tree_->cleanTree(this,-COIN_DBL_MAX,bestPossibleObjective_) ;
+      delete nextRowCut_;
+      if (stoppedOnGap_)
+        { messageHandler()->message(CBC_GAP,messages())
+          << bestObjective_-bestPossibleObjective_
+          << dblParam_[CbcAllowableGap]
+          << dblParam_[CbcAllowableFractionGap]*100.0
+          << CoinMessageEol ;
+        secondaryStatus_ = 2;
+        status_ = 0 ; }
+        else
+          if (isNodeLimitReached())
+            { handler_->message(CBC_MAXNODES,messages_) << CoinMessageEol ;
+            secondaryStatus_ = 3;
+            status_ = 1 ; }
+          else
+        if (totalTime >= dblParam_[CbcMaximumSeconds])
+          { handler_->message(CBC_MAXTIME,messages_) << CoinMessageEol ; 
+          secondaryStatus_ = 4;
+          status_ = 1 ; }
+        else
+          if (eventHappened_)
+            { handler_->message(CBC_EVENT,messages_) << CoinMessageEol ; 
+            secondaryStatus_ = 5;
+            status_ = 5 ; }
+          else
+            { handler_->message(CBC_MAXSOLS,messages_) << CoinMessageEol ;
+            secondaryStatus_ = 6;
+            status_ = 1 ; }
+    }
 /*
   That's it, we've exhausted the search tree, or broken out of the loop because
@@ -2464 +2780 @@
   subTreeModel_(NULL),
   numberStoppedSubTrees_(0),
+  mutex_(NULL),
   presolve_(0),
   numberStrong_(5),
@@ -2511 +2828 @@
   searchStrategy_(-1),
   numberStrongIterations_(0),
-  resolveAfterTakeOffCuts_(true)
+  resolveAfterTakeOffCuts_(true),
+#if NEW_UPDATE_OBJECT>1
+  numberUpdateItems_(0),
+  maximumNumberUpdateItems_(0),
+  updateItems_(NULL),
+#endif
+  numberThreads_(0),
+  threadMode_(0)
 {
   memset(intParam_,0,sizeof(intParam_));
@@ -2594 +2918 @@
   subTreeModel_(NULL),
   numberStoppedSubTrees_(0),
+  mutex_(NULL),
   presolve_(0),
   numberStrong_(5),
@@ -2641 +2966 @@
   searchStrategy_(-1),
   numberStrongIterations_(0),
-  resolveAfterTakeOffCuts_(true)
+  resolveAfterTakeOffCuts_(true),
+#if NEW_UPDATE_OBJECT>1
+  numberUpdateItems_(0),
+  maximumNumberUpdateItems_(0),
+  updateItems_(NULL),
+#endif
+  numberThreads_(0),
+  threadMode_(0)
 {
   memset(intParam_,0,sizeof(intParam_));
@@ -2780 +3112 @@
       numberIntegers_++;
   }
+  delete [] integerVariable_;
   if (numberIntegers_) {
-    delete [] integerVariable_;
     integerVariable_ = new int [numberIntegers_];
     numberIntegers_=0;
@@ -2819 +3151 @@
   subTreeModel_(rhs.subTreeModel_),
   numberStoppedSubTrees_(rhs.numberStoppedSubTrees_),
+  mutex_(NULL),
   presolve_(rhs.presolve_),
   numberStrong_(rhs.numberStrong_),
@@ -2855 +3188 @@
   searchStrategy_(rhs.searchStrategy_),
   numberStrongIterations_(rhs.numberStrongIterations_),
-  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_)
+  resolveAfterTakeOffCuts_(rhs.resolveAfterTakeOffCuts_),
+#if NEW_UPDATE_OBJECT>1
+  numberUpdateItems_(rhs.numberUpdateItems_),
+  maximumNumberUpdateItems_(rhs.maximumNumberUpdateItems_),
+  updateItems_(NULL),
+#endif
+  numberThreads_(rhs.numberThreads_),
+  threadMode_(rhs.threadMode_)
 {
   memcpy(intParam_,rhs.intParam_,sizeof(intParam_));
@@ -2907 +3247 @@
       object_ = new OsiObject * [numberObjects_];
       int i;
-      for (i=0;i<numberObjects_;i++) 
+      for (i=0;i<numberObjects_;i++) {
         object_[i]=(rhs.object_[i])->clone();
+        CbcObject * obj = dynamic_cast <CbcObject *>(object_[i]) ;
+        // Could be OsiObjects
+        if (obj)
+          obj->setModel(this);
+      }
     } else {
       object_=NULL;
@@ -2932 +3277 @@
     originalColumns_=NULL;
   }
+#if NEW_UPDATE_OBJECT>1
+  if (maximumNumberUpdateItems_) {
+    updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
+    for (int i=0;i<maximumNumberUpdateItems_;i++)
+      updateItems_[i] = rhs.updateItems_[i];
+  }
+#endif
   if (maximumWhich_&&rhs.whichGenerator_)
     whichGenerator_ = CoinCopyOfArray(rhs.whichGenerator_,maximumWhich_);
@@ -3105 +3457 @@
     subTreeModel_ = rhs.subTreeModel_;
     numberStoppedSubTrees_ = rhs.numberStoppedSubTrees_;
+    mutex_ = NULL;
     presolve_ = rhs.presolve_;
     numberStrong_ = rhs.numberStrong_;
@@ -3156 +3509 @@
     numberNewCuts_ = rhs.numberNewCuts_;
     resolveAfterTakeOffCuts_=rhs.resolveAfterTakeOffCuts_;
+#if NEW_UPDATE_OBJECT>1
+    numberUpdateItems_ = rhs.numberUpdateItems_;
+    maximumNumberUpdateItems_ = rhs.maximumNumberUpdateItems_;
+    delete [] updateItems_;
+    if (maximumNumberUpdateItems_) {
+      updateItems_ = new CbcObjectUpdateData [maximumNumberUpdateItems_];
+      for (i=0;i<maximumNumberUpdateItems_;i++)
+        updateItems_[i] = rhs.updateItems_[i];
+    } else {
+      updateItems_ = NULL;
+    }
+#endif
+    numberThreads_ = rhs.numberThreads_;
+    threadMode_ = rhs.threadMode_;
     sizeMiniTree_ = rhs.sizeMiniTree_;
     searchStrategy_ = rhs.searchStrategy_;
@@ -3329 +3696 @@
   originalColumns_=NULL;
   delete strategy_;
+#if NEW_UPDATE_OBJECT>1
+  delete [] updateItems_;
+  updateItems_=NULL;
+  numberUpdateItems_=0;
+  maximumNumberUpdateItems_=0;
+#endif
   gutsOfDestructor2();
 }
@@ -3686 +4059 @@
 */
   currentNumberCuts_=currentNumberCuts;
-  if (currentNumberCuts >= maximumNumberCuts_) {
+  if (currentNumberCuts > maximumNumberCuts_) {
     maximumNumberCuts_ = currentNumberCuts;
     delete [] addedCuts_;
@@ -3784 +4157 @@
   reconstructed basis in the solver.
 */
-  if (node->objectiveValue() < cutoff)
+  if (node->objectiveValue() < cutoff||numberThreads_)
   { 
 #   ifdef CBC_CHECK_BASIS
@@ -3864 +4237 @@
   set. Adjust the cut reference counts to reflect that we no longer need to
   explore the remaining branch arms, hence they will no longer reference any
-  cuts. Cuts whose reference count falls to zero are deleted.
+  cuts. Cuts whose reference count falls to zero are deleted.  
 */
   else
   { int i;
-    int numberLeft = nodeInfo->numberBranchesLeft();
-    for (i = 0 ; i < currentNumberCuts ; i++)
-    { if (addedCuts_[i])
-      { if (!addedCuts_[i]->decrement(numberLeft))
-        { delete addedCuts_[i];
-          addedCuts_[i] = NULL; } } }
+    if (currentNumberCuts) {
+      lockThread();
+      int numberLeft = nodeInfo->numberBranchesLeft();
+      for (i = 0 ; i < currentNumberCuts ; i++)
+        { if (addedCuts_[i])
+          { if (!addedCuts_[i]->decrement(numberLeft))
+            { delete addedCuts_[i];
+            addedCuts_[i] = NULL; } } }
+      unlockThread();
+    }
     return 1 ; }
 }
@@ -3988 +4365 @@
     saveClpOptions = clpSolver->specialOptions();
 # endif
+#ifdef CBC_THREAD
+  CbcModel ** threadModel = NULL;
+  pthread_t * threadId = NULL;
+  pthread_cond_t condition_main;
+  pthread_mutex_t condition_mutex;
+  pthread_mutex_t * mutex2 = NULL;
+  pthread_cond_t * condition2 = NULL;
+  threadStruct * threadInfo = NULL;
+  void * saveMutex = NULL;
+  if (numberThreads_&&(threadMode_&2)!=0&&!numberNodes_) {
+    threadId = new pthread_t [numberThreads_];
+    pthread_cond_init(&condition_main,NULL);
+    pthread_mutex_init(&condition_mutex,NULL);
+    threadModel = new CbcModel * [numberThreads_];
+    threadInfo = new threadStruct [numberThreads_+1];
+    mutex2 = new pthread_mutex_t [numberThreads_];
+    condition2 = new pthread_cond_t [numberThreads_];
+    saveMutex = mutex_;
+    for (int i=0;i<numberThreads_;i++) {
+      pthread_mutex_init(mutex2+i,NULL);
+      pthread_cond_init(condition2+i,NULL);
+      threadId[i]=0;
+      threadModel[i]=new CbcModel;
+      threadModel[i]->generator_ = new CbcCutGenerator * [1];
+      delete threadModel[i]->solver_;
+      threadModel[i]->solver_=NULL;
+      threadModel[i]->numberThreads_=numberThreads_;
+      mutex_ = (void *) (threadInfo+i);
+      threadInfo[i].thisModel=(CbcModel *) threadModel[i];
+      threadInfo[i].baseModel=this;
+      threadInfo[i].threadIdOfBase=pthread_self();
+      threadInfo[i].mutex2=mutex2+i;
+      threadInfo[i].condition2=condition2+i;
+      threadInfo[i].returnCode=-1;
+      pthread_create(threadId+i,NULL,doCutsThread,threadInfo+i);
+    }
+    // Do a partial one for base model
+    threadInfo[numberThreads_].baseModel=this;
+    mutex_ = (void *) (threadInfo+numberThreads_);
+    threadInfo[numberThreads_].condition2=&condition_main;
+    threadInfo[numberThreads_].mutex2=&condition_mutex;
+  }
+#endif
   bool feasible = true ;
   int lastNumberCuts = 0 ;
@@ -4026 +4446 @@
   //printf("zero iterations on first solve of branch\n");
 #endif
-  if (node&&!node->nodeInfo()->numberBranchesLeft())
+  if (node&&node->nodeInfo()&&!node->nodeInfo()->numberBranchesLeft())
     node->nodeInfo()->allBranchesGone(); // can clean up
   feasible = returnCode  != 0 ;
@@ -4035 +4455 @@
   }
   
+#if NEW_UPDATE_OBJECT==0
   // Update branching information if wanted
   if(node &&branchingMethod_)
     branchingMethod_->updateInformation(solver_,node);
+#elif NEW_UPDATE_OBJECT<2
+  // Update branching information if wanted
+  if(node &&branchingMethod_) {
+    OsiBranchingObject * bobj = node->modifiableBranchingObject();
+    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
+    if (cbcobj) {
+      CbcObject * object = cbcobj->object();
+      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
+      object->updateInformation(update);
+    } else {
+      branchingMethod_->updateInformation(solver_,node);
+    }
+  }
+#else
+  // Update branching information if wanted
+  if(node &&branchingMethod_) {
+    OsiBranchingObject * bobj = node->modifiableBranchingObject();
+    CbcBranchingObject * cbcobj = dynamic_cast<CbcBranchingObject *> (bobj);
+    if (cbcobj) {
+      CbcObject * object = cbcobj->object();
+      CbcObjectUpdateData update = object->createUpdateInformation(solver_,node,cbcobj);
+      // have to compute object number as not saved
+      CbcSimpleInteger * simpleObject =
+          dynamic_cast <CbcSimpleInteger *>(object) ;
+      int iObject;
+      int iColumn = simpleObject->columnNumber();
+      for (iObject = 0 ; iObject < numberObjects_ ; iObject++) {
+        simpleObject =
+          dynamic_cast <CbcSimpleInteger *>(object_[iObject]) ;
+        if (simpleObject->columnNumber()==iColumn)
+          break;
+      }
+      assert (iObject<numberObjects_);
+      update.objectNumber_ = iObject;
+      addUpdateInformation(update);
+    } else {
+      branchingMethod_->updateInformation(solver_,node);
+    }
+  }
+#endif
 
 #ifdef CBC_DEBUG
@@ -4078 +4539 @@
   sumChangeObjective1_ += solver_->getObjValue()*solver_->getObjSense()
     - objectiveValue ;
-  if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
+  if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
     numberTries=0; // exit
   //if ((numberNodes_%100)==0)
@@ -4244 +4705 @@
       whichGenerator_[numberViolated++]=-2;
     }
-    double * newSolution = new double [numberColumns] ;
-    double heuristicValue = getCutoff() ;
-    int found = -1; // no solution found
 
     // reset probing info
     if (probingInfo_)
       probingInfo_->initializeFixing();
-    for (int i = 0;i<numberCutGenerators_+numberHeuristics_;i++) {
-      int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
-      int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
-      if (i<numberCutGenerators_) {
-        bool generate = generator_[i]->normal();
-        // skip if not optimal and should be (maybe a cut generator has fixed variables)
-        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
-          generate=false;
-        if (generator_[i]->switchedOff())
-          generate=false;;
+    int i;
+    if ((threadMode_&2)==0||numberNodes_) {
+      for (i = 0;i<numberCutGenerators_;i++) {
+        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
+        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
+        bool generate = generator_[i]->normal();
+        // skip if not optimal and should be (maybe a cut generator has fixed variables)
+        if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
+          generate=false;
+        if (generator_[i]->switchedOff())
+          generate=false;;
         if (generate) {
           bool mustResolve = 
-            generator_[i]->generateCuts(theseCuts,fullScan,node) ;
+            generator_[i]->generateCuts(theseCuts,fullScan,solver_,node) ;
           int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
-          if(numberRowCutsBefore < numberRowCutsAfter &&
-             generator_[i]->mustCallAgain())
-            keepGoing=true; // say must go round
+          if(numberRowCutsBefore < numberRowCutsAfter &&
+             generator_[i]->mustCallAgain())
+            keepGoing=true; // say must go round
           // Check last cut to see if infeasible
-          if(numberRowCutsBefore < numberRowCutsAfter) {
+          if(numberRowCutsBefore < numberRowCutsAfter) {
             const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
             if (thisCut->lb()>thisCut->ub()) {
@@ -4277 +4736 @@
           }
 #ifdef CBC_DEBUG
-          {
-            int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
-            int k ;
-            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
-              OsiRowCut thisCut = theseCuts.rowCut(k) ;
-              /* check size of elements.
-                 We can allow smaller but this helps debug generators as it
-                 is unsafe to have small elements */
-              int n=thisCut.row().getNumElements();
-              const int * column = thisCut.row().getIndices();
-              const double * element = thisCut.row().getElements();
-              //assert (n);
-              for (int i=0;i<n;i++) {
-                double value = element[i];
-                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
-              }
-            }
-          }
+          {
+            int k ;
+            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
+              OsiRowCut thisCut = theseCuts.rowCut(k) ;
+              /* check size of elements.
+                 We can allow smaller but this helps debug generators as it
+                 is unsafe to have small elements */
+              int n=thisCut.row().getNumElements();
+              const int * column = thisCut.row().getIndices();
+              const double * element = thisCut.row().getElements();
+              //assert (n);
+              for (int i=0;i<n;i++) {
+                double value = element[i];
+                assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
+              }
+            }
+          }
 #endif
           if (mustResolve) {
-            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
-            feasible = returnCode  != 0 ;
-            if (returncode<0)
-              numberTries=0;
-            if ((specialOptions_&1)!=0) {
-              debugger = solver_->getRowCutDebugger() ;
-              if (debugger) 
-                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
-              else
-                onOptimalPath=false;
-              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
-                assert(feasible) ;
-            }
+            int returncode = resolve(node ? node->nodeInfo() : NULL,2);
+            feasible = returnCode  != 0 ;
+            if (returncode<0)
+              numberTries=0;
+            if ((specialOptions_&1)!=0) {
+              debugger = solver_->getRowCutDebugger() ;
+              if (debugger) 
+                onOptimalPath = (debugger->onOptimalPath(*solver_)) ;
+              else
+                onOptimalPath=false;
+              if (onOptimalPath && !solver_->isDualObjectiveLimitReached())
+                assert(feasible) ;
+            }
             if (!feasible)
               break ;
           }
         }
+        int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
+        int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
+        
+        if ((specialOptions_&1)!=0) {
+          if (onOptimalPath) {
+            int k ;
+            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
+              OsiRowCut thisCut = theseCuts.rowCut(k) ;
+              if(debugger->invalidCut(thisCut)) {
+                solver_->writeMps("badCut");
+#ifdef NDEBUG
+                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
+                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
+                abort();
+#endif
+              }
+              assert(!debugger->invalidCut(thisCut)) ;
+            }
+          }
+        }
+/*
+  The cut generator has done its thing, and maybe it generated some
+  cuts.  Do a bit of bookkeeping: load
+  whichGenerator[i] with the index of the generator responsible for a cut,
+  and place cuts flagged as global in the global cut pool for the model.
+
+  lastNumberCuts is the sum of cuts added in previous iterations; it's the
+  offset to the proper starting position in whichGenerator.
+*/
+        int numberBefore =
+          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
+        int numberAfter =
+          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
+        // possibly extend whichGenerator
+        resizeWhichGenerator(numberBefore, numberAfter);
+        int j ;
+        if (fullScan) {
+          // counts
+          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
+        }
+        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
+        
+        bool dodgyCuts=false;
+        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
+          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
+          if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) {
+            dodgyCuts=true;
+            break;
+          }
+          whichGenerator_[numberBefore++] = i ;
+          if (thisCut->lb()>thisCut->ub())
+            violated=-2; // sub-problem is infeasible
+          if (thisCut->globallyValid()) {
+            // add to global list
+            OsiRowCut newCut(*thisCut);
+            newCut.setGloballyValid(2);
+            newCut.mutableRow().setTestForDuplicateIndex(false);
+            globalCuts_.insert(newCut) ;
+          }
+        }
+        if (dodgyCuts) {
+          for (int k=numberRowCutsAfter-1;k>=j;k--) {
+            const OsiRowCut * thisCut = theseCuts.rowCutPtr(k) ;
+            if (thisCut->lb()>thisCut->ub())
+              violated=-2; // sub-problem is infeasible
+            if (thisCut->lb()>1.0e10||thisCut->ub()<-1.0e10) 
+              theseCuts.eraseRowCut(k);
+          }
+          numberRowCutsAfter = theseCuts.sizeRowCuts() ;
+          for (;j<numberRowCutsAfter;j++) {
+            const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
+            whichGenerator_[numberBefore++] = i ;
+            if (thisCut->globallyValid()) {
+              // add to global list
+              OsiRowCut newCut(*thisCut);
+              newCut.setGloballyValid(2);
+              newCut.mutableRow().setTestForDuplicateIndex(false);
+              globalCuts_.insert(newCut) ;
+            }
+          }
+        }
+        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
+          whichGenerator_[numberBefore++] = i ;
+          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
+          if (thisCut->globallyValid()) {
+            // add to global list
+            OsiColCut newCut(*thisCut);
+            newCut.setGloballyValid(2);
+            globalCuts_.insert(newCut) ;
+          }
+        }
+      }
+      // Add in any violated saved cuts
+      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
+        int numberOld = theseCuts.sizeRowCuts();
+        int numberCuts = slackCuts.sizeRowCuts() ;
+        int i;
+        // possibly extend whichGenerator
+        resizeWhichGenerator(numberOld, numberOld+numberCuts);
+        for ( i = 0;i<numberCuts;i++) {
+          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
+          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
+            if (messageHandler()->logLevel()>2)
+              printf("Old cut added - violation %g\n",
+                     thisCut->violated(cbcColSolution_)) ;
+            whichGenerator_[numberOld++]=-1;
+            theseCuts.insert(*thisCut) ;
+          }
+        }
+      }
+    } else {
+      // do cuts independently
+      OsiCuts * eachCuts = new OsiCuts [numberCutGenerators_];;
+#ifdef CBC_THREAD
+      if (!threadModel) {
+#endif
+        // generate cuts
+        for (i = 0;i<numberCutGenerators_;i++) {
+          bool generate = generator_[i]->normal();
+          // skip if not optimal and should be (maybe a cut generator has fixed variables)
+          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
+            generate=false;
+          if (generator_[i]->switchedOff())
+            generate=false;;
+          if (generate) 
+            generator_[i]->generateCuts(eachCuts[i],fullScan,solver_,node) ;
+        }
+#ifdef CBC_THREAD
       } else {
+        for (i=0;i<numberThreads_;i++) {
+          // set solver here after cloning
+          threadModel[i]->solver_=solver_->clone();
+          threadModel[i]->numberNodes_ = (fullScan) ? 1 : 0;
+        }
+        // generate cuts
+        for (i = 0;i<numberCutGenerators_;i++) {
+          bool generate = generator_[i]->normal();
+          // skip if not optimal and should be (maybe a cut generator has fixed variables)
+          if (generator_[i]->needsOptimalBasis()&&!solver_->basisIsAvailable())
+            generate=false;
+          if (generator_[i]->switchedOff())
+            generate=false;;
+          if (generate) { 
+            bool finished=false;
+            int iThread=-1;
+            // see if any available
+            for (iThread=0;iThread<numberThreads_;iThread++) {
+              if (threadInfo[iThread].returnCode) {
+                finished=true;
+                break;
+              } else if (threadInfo[iThread].returnCode==0) {
+                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+              }
+            }
+            while (!finished) {
+              pthread_mutex_lock(&condition_mutex);
+              struct timespec absTime;
+              clock_gettime(CLOCK_REALTIME,&absTime);
+              absTime.tv_nsec += 1000000; // millisecond
+              if (absTime.tv_nsec>=1000000000) {
+                absTime.tv_nsec -= 1000000000;
+                absTime.tv_sec++;
+              }
+              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
+              pthread_mutex_unlock(&condition_mutex);
+              for (iThread=0;iThread<numberThreads_;iThread++) {
+                if (threadInfo[iThread].returnCode>0) {
+                  finished=true;
+                  break;
+                } else if (threadInfo[iThread].returnCode==0) {
+                  pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+                }
+              }
+            }
+            assert (iThread<numberThreads_);
+            assert (threadInfo[iThread].returnCode);
+            threadModel[iThread]->generator_[0]=generator_[i];
+            threadModel[iThread]->object_ = (OsiObject **) (eachCuts+i);
+            // allow to start
+            threadInfo[iThread].returnCode=0;
+            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+          }
+        }
+        // wait
+        for (int iThread=0;iThread<numberThreads_;iThread++) {
+          if (threadInfo[iThread].returnCode==0) {
+            bool finished=false;
+            pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+            while (!finished) {
+              pthread_mutex_lock(&condition_mutex);
+              struct timespec absTime;
+              clock_gettime(CLOCK_REALTIME,&absTime);
+              absTime.tv_nsec += 1000000; // millisecond
+              if (absTime.tv_nsec>=1000000000) {
+                absTime.tv_nsec -= 1000000000;
+                absTime.tv_sec++;
+              }
+              pthread_cond_timedwait(&condition_main,&condition_mutex,&absTime);
+              pthread_mutex_unlock(&condition_mutex);
+              if (threadInfo[iThread].returnCode>0) {
+                finished=true;
+                break;
+              } else if (threadInfo[iThread].returnCode==0) {
+                pthread_cond_signal(threadInfo[iThread].condition2); // unlock
+              }
+            }
+          }
+          assert (threadInfo[iThread].returnCode);
+          // say available
+          threadInfo[iThread].returnCode=-1;
+          delete threadModel[iThread]->solver_;
+          threadModel[iThread]->solver_=NULL;
+        }
+      }
+#endif
+      // Now put together
+      for (i = 0;i<numberCutGenerators_;i++) {
+        // add column cuts
+        int numberColumnCutsBefore = theseCuts.sizeColCuts() ;
+        int numberColumnCuts = eachCuts[i].sizeColCuts();
+        int numberColumnCutsAfter = numberColumnCutsBefore
+          + numberColumnCuts;
+        int j;
+        for (j=0;j<numberColumnCuts;j++) {
+          theseCuts.insert(eachCuts[i].colCut(j));
+        }
+        int numberRowCutsBefore = theseCuts.sizeRowCuts() ;
+        int numberRowCuts = eachCuts[i].sizeRowCuts();
+        int numberRowCutsAfter = numberRowCutsBefore
+          + numberRowCuts;
+        if (numberRowCuts) {
+          for (j=0;j<numberRowCuts;j++) {
+            const OsiRowCut * thisCut = eachCuts[i].rowCutPtr(j) ;
+            if (thisCut->lb()<=1.0e10&&thisCut->ub()>=-1.0e10) 
+              theseCuts.insert(eachCuts[i].rowCut(j));
+          }
+          if (generator_[i]->mustCallAgain())
+            keepGoing=true; // say must go round
+          // Check last cut to see if infeasible
+          const OsiRowCut * thisCut = theseCuts.rowCutPtr(numberRowCutsAfter-1) ;
+          if (thisCut->lb()>thisCut->ub()) {
+            feasible = false; // sub-problem is infeasible
+            break;
+          }
+        }
+#ifdef CBC_DEBUG
+        {
+          int k ;
+          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
+            OsiRowCut thisCut = theseCuts.rowCut(k) ;
+            /* check size of elements.
+               We can allow smaller but this helps debug generators as it
+               is unsafe to have small elements */
+            int n=thisCut.row().getNumElements();
+            const int * column = thisCut.row().getIndices();
+            const double * element = thisCut.row().getElements();
+            //assert (n);
+            for (int i=0;i<n;i++) {
+              double value = element[i];
+              assert(fabs(value)>1.0e-12&&fabs(value)<1.0e20);
+            }
+          }
+        }
+#endif
+        if ((specialOptions_&1)!=0) {
+          if (onOptimalPath) {
+            int k ;
+            for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
+              OsiRowCut thisCut = theseCuts.rowCut(k) ;
+              if(debugger->invalidCut(thisCut)) {
+                solver_->writeMps("badCut");
+#ifdef NDEBUG
+                printf("Cut generator %d (%s) produced invalid cut (%dth in this go)\n",
+                       i,generator_[i]->cutGeneratorName(),k-numberRowCutsBefore);
+                abort();
+#endif
+              }
+              assert(!debugger->invalidCut(thisCut)) ;
+            }
+          }
+        }
+/*
+  The cut generator has done its thing, and maybe it generated some
+  cuts.  Do a bit of bookkeeping: load
+  whichGenerator[i] with the index of the generator responsible for a cut,
+  and place cuts flagged as global in the global cut pool for the model.
+
+  lastNumberCuts is the sum of cuts added in previous iterations; it's the
+  offset to the proper starting position in whichGenerator.
+*/
+        int numberBefore =
+          numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
+        int numberAfter =
+          numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
+        // possibly extend whichGenerator
+        resizeWhichGenerator(numberBefore, numberAfter);
+        if (fullScan) {
+          // counts
+          countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
+        }
+        countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
+        
+        for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
+          whichGenerator_[numberBefore++] = i ;
+          const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
+          if (thisCut->lb()>thisCut->ub())
+            violated=-2; // sub-problem is infeasible
+          if (thisCut->globallyValid()) {
+            // add to global list
+            OsiRowCut newCut(*thisCut);
+            newCut.setGloballyValid(2);
+            newCut.mutableRow().setTestForDuplicateIndex(false);
+            globalCuts_.insert(newCut) ;
+          }
+        }
+        for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
+          whichGenerator_[numberBefore++] = i ;
+          const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
+          if (thisCut->globallyValid()) {
+            // add to global list
+            OsiColCut newCut(*thisCut);
+            newCut.setGloballyValid(2);
+            globalCuts_.insert(newCut) ;
+          }
+        }
+      }
+      // Add in any violated saved cuts
+      if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
+        int numberOld = theseCuts.sizeRowCuts();
+        int numberCuts = slackCuts.sizeRowCuts() ;
+        int i;
+        // possibly extend whichGenerator
+        resizeWhichGenerator(numberOld, numberOld+numberCuts);
+        for ( i = 0;i<numberCuts;i++) {
+          const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
+          if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
+            if (messageHandler()->logLevel()>2)
+              printf("Old cut added - violation %g\n",
+                     thisCut->violated(cbcColSolution_)) ;
+            whichGenerator_[numberOld++]=-1;
+            theseCuts.insert(*thisCut) ;
+          }
+        }
+      }
+      delete [] eachCuts;
+    }
+    //if (!feasible)
+    //break;
+/*
+  End of the loop to exercise each generator - try heuristics
+  - unless at root node and first pass
+*/
+    if (numberNodes_||currentPassNumber_!=1) {
+      double * newSolution = new double [numberColumns] ;
+      double heuristicValue = getCutoff() ;
+      int found = -1; // no solution found
+      for (i = 0;i<numberHeuristics_;i++) {
         // see if heuristic will do anything
         double saveValue = heuristicValue ;
         int ifSol = 
-          heuristic_[i-numberCutGenerators_]->solution(heuristicValue,
-                                                       newSolution,
-                                                       theseCuts) ;
+          heuristic_[i]->solution(heuristicValue,
+                                  newSolution,
+                                  theseCuts) ;
         if (ifSol>0) {
           // better solution found
-          found = i-numberCutGenerators_ ;
+          found = i ;
           incrementUsed(newSolution);
         } else if (ifSol<0) {
@@ -4329 +5143 @@
         }
       }
-      int numberRowCutsAfter = theseCuts.sizeRowCuts() ;
-      int numberColumnCutsAfter = theseCuts.sizeColCuts() ;
-
-      if ((specialOptions_&1)!=0) {
-        if (onOptimalPath) {
-          int k ;
-          for (k = numberRowCutsBefore;k<numberRowCutsAfter;k++) {
-            OsiRowCut thisCut = theseCuts.rowCut(k) ;
-            if(debugger->invalidCut(thisCut)) {
-              solver_->writeMps("badCut");
-            }
-            assert(!debugger->invalidCut(thisCut)) ;
-          }
-        }
-      }
-/*
-  The cut generator/heuristic has done its thing, and maybe it generated some
-  cuts and/or a new solution.  Do a bit of bookkeeping: load
-  whichGenerator[i] with the index of the generator responsible for a cut,
-  and place cuts flagged as global in the global cut pool for the model.
-
-  lastNumberCuts is the sum of cuts added in previous iterations; it's the
-  offset to the proper starting position in whichGenerator.
-*/
-      int numberBefore =
-            numberRowCutsBefore+numberColumnCutsBefore+lastNumberCuts ;
-      int numberAfter =
-            numberRowCutsAfter+numberColumnCutsAfter+lastNumberCuts ;
-      // possibly extend whichGenerator
-      resizeWhichGenerator(numberBefore, numberAfter);
-      int j ;
-      if (fullScan) {
-        // counts
-        countColumnCuts[i] += numberColumnCutsAfter-numberColumnCutsBefore ;
-      }
-      countRowCuts[i] += numberRowCutsAfter-numberRowCutsBefore ;
-      
-      for (j = numberRowCutsBefore;j<numberRowCutsAfter;j++) {
-        whichGenerator_[numberBefore++] = i ;
-        const OsiRowCut * thisCut = theseCuts.rowCutPtr(j) ;
-        if (thisCut->lb()>thisCut->ub())
-          violated=-2; // sub-problem is infeasible
-        if (thisCut->globallyValid()) {
-          // add to global list
-          OsiRowCut newCut(*thisCut);
-          newCut.setGloballyValid(2);
-          newCut.mutableRow().setTestForDuplicateIndex(false);
-          globalCuts_.insert(newCut) ;
-        }
-      }
-      for (j = numberColumnCutsBefore;j<numberColumnCutsAfter;j++) {
-        whichGenerator_[numberBefore++] = i ;
-        const OsiColCut * thisCut = theseCuts.colCutPtr(j) ;
-        if (thisCut->globallyValid()) {
-          // add to global list
-          OsiColCut newCut(*thisCut);
-          newCut.setGloballyValid(2);
-          globalCuts_.insert(newCut) ;
-        }
-      }
-    }
-    // If at root node and first pass do heuristics without cuts
-    if (!numberNodes_&&currentPassNumber_==1) {
-      // Save number solutions
-      int saveNumberSolutions = numberSolutions_;
-      int saveNumberHeuristicSolutions = numberHeuristicSolutions_;
-      // make sure bounds tight
-      { 
-        for (int i = 0;i < numberObjects_;i++)
-          object_[i]->resetBounds(solver_) ;
-      }
-      for (int i = 0;i<numberHeuristics_;i++) {
-        // see if heuristic will do anything
-        double saveValue = heuristicValue ;
-        int ifSol = heuristic_[i]->solution(heuristicValue,
-                                            newSolution);
-        if (ifSol>0) {
-          // better solution found
-          found = i ;
-          incrementUsed(newSolution);
-          // increment number of solutions so other heuristics can test
-          numberSolutions_++;
-          numberHeuristicSolutions_++;
-        } else {
-          heuristicValue = saveValue ;
-        }
-      }
-      // Restore number solutions
-      numberSolutions_ = saveNumberSolutions;
-      numberHeuristicSolutions_ = saveNumberHeuristicSolutions;
-    }
-    // Add in any violated saved cuts
-    if (!theseCuts.sizeRowCuts()&&!theseCuts.sizeColCuts()) {
-      int numberOld = theseCuts.sizeRowCuts();
-      int numberCuts = slackCuts.sizeRowCuts() ;
-      int i;
-      // possibly extend whichGenerator
-      resizeWhichGenerator(numberOld, numberOld+numberCuts);
-      for ( i = 0;i<numberCuts;i++) {
-        const OsiRowCut * thisCut = slackCuts.rowCutPtr(i) ;
-        if (thisCut->violated(cbcColSolution_)>100.0*primalTolerance) {
-          if (messageHandler()->logLevel()>2)
-            printf("Old cut added - violation %g\n",
-                   thisCut->violated(cbcColSolution_)) ;
-          whichGenerator_[numberOld++]=-1;
-          theseCuts.insert(*thisCut) ;
-        }
-      }
-    }
-/*
-  End of the loop to exercise each generator/heuristic.
-
+/*
   Did any of the heuristics turn up a new solution? Record it before we free
   the vector.
 */
-    if (found >= 0) { 
-      phase_=4;
-      incrementUsed(newSolution);
-      setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
-      lastHeuristic_ = heuristic_[found];
-      CbcTreeLocal * tree 
+      if (found >= 0) { 
+        phase_=4;
+        incrementUsed(newSolution);
+        lastHeuristic_ = heuristic_[found];
+        setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
+        CbcTreeLocal * tree 
           = dynamic_cast<CbcTreeLocal *> (tree_);
-      if (tree)
-        tree->passInSolution(bestSolution_,heuristicValue);
-    }
-    delete [] newSolution ;
+        if (tree)
+          tree->passInSolution(bestSolution_,heuristicValue);
+      }
+      delete [] newSolution ;
+    }
 
 #if 0
@@ -4616 +5320 @@
       }
       feasible = resolve(node ? node->nodeInfo() : NULL,2) ;
-      if ( CoinCpuTime()-dblParam_[CbcStartSeconds] > dblParam_[CbcMaximumSeconds] )
+      if ( getCurrentSeconds() > dblParam_[CbcMaximumSeconds] )
         numberTries=0; // exit
 #     ifdef CBC_DEBUG
@@ -4689 +5393 @@
     if (!feasible)
     { int i ;
-      for (i = 0;i<currentNumberCuts_;i++) {
-        // take off node
-        if (addedCuts_[i]) {
-          if (!addedCuts_[i]->decrement())
-            delete addedCuts_[i] ;
-          addedCuts_[i] = NULL ;
+      if (currentNumberCuts_) {
+        lockThread();
+        for (i = 0;i<currentNumberCuts_;i++) {
+          // take off node
+          if (addedCuts_[i]) {
+            if (!addedCuts_[i]->decrement())
+              delete addedCuts_[i] ;
+            addedCuts_[i] = NULL ;
+          }
         }
       }
@@ -4830 +5537 @@
       phase_=4;
       incrementUsed(newSolution);
+      lastHeuristic_ = heuristic_[found];
       setBestSolution(CBC_ROUNDING,heuristicValue,newSolution) ;
-      lastHeuristic_ = heuristic_[found];
     }
     delete [] newSolution ;
@@ -5165 +5872 @@
     clpSolver->setSpecialOptions(saveClpOptions);
 # endif
+#ifdef CBC_THREAD
+  if (threadModel) {
+    // stop threads
+    int i;
+    for (i=0;i<numberThreads_;i++) {
+      while (threadInfo[i].returnCode==0) {
+        pthread_cond_signal(threadInfo[i].condition2); // unlock
+        pthread_mutex_lock(&condition_mutex);
+        struct timespec absTime;