Add parallel lock. Fix cut pool sepa

Author: Opt-Mucca

highs/mip/HighsCutPool.cpp

@@ -148,18 +148,21 @@ double HighsCutPool::getParallelism(HighsInt row1, HighsInt row2,
 }
 
 void HighsCutPool::lpCutRemoved(HighsInt cut, bool thread_safe) {
-  numLps_[cut].fetch_add(-1, std::memory_order_relaxed);
+  HighsInt numLps = numLps_[cut].fetch_add(-1, std::memory_order_relaxed);
   if (thread_safe) return;
   if (matrix_.columnsLinked(cut)) {
-    propRows.erase(std::make_pair(-1, cut));
+    propRows.erase(std::make_pair(ages_[cut], cut));
     propRows.emplace(1, cut);
   }
   ages_[cut] = 1;
   --numLpCuts;
   ++ageDistribution[1];
+  if (numLps == 0) {
+    numLps_[cut].fetch_add(-1, std::memory_order_relaxed);
+  }
 }
 
-void HighsCutPool::performAging(const bool parallel_sepa) {
+void HighsCutPool::performAging() {
   HighsInt cutIndexEnd = matrix_.getNumRows();
 
   HighsInt agelim = agelim_;
@@ -170,19 +173,17 @@ void HighsCutPool::performAging(const bool parallel_sepa) {
   }
 
   for (HighsInt i = 0; i != cutIndexEnd; ++i) {
-    if (numLps_[i] == 0 && ages_[i] >= 0 && parallel_sepa) {
-      ageDistribution[ages_[i]] -= 1;
-      lpCutRemoved(i);
-      numLps_[i] = -1;
-    } else if (numLps_[i] > 0 && ages_[i] >= 0 && parallel_sepa) {
-      // Age and propRows were not updated in the multi-thread case
+    if (numLps_[i] > 0 && ages_[i] >= 0) {
+      --ageDistribution[ages_[i]];
       if (matrix_.columnsLinked(i)) {
         propRows.erase(std::make_pair(ages_[i], i));
         propRows.emplace(-1, i);
       }
-      --ageDistribution[ages_[i]];
       ages_[i] = -1;
     }
+    if (numLps_[i] == 0) {
+      lpCutRemoved(i);
+    }
     if (ages_[i] < 0) continue;
 
     bool isPropagated = matrix_.columnsLinked(i);
@@ -369,7 +370,9 @@ void HighsCutPool::separate(const std::vector<double>& sol, HighsDomain& domain,
 
   efficacious_cuts.resize(numefficacious);
 
-  HighsInt selectednnz = cutset.ARindex_.size();
+  HighsInt orignumcuts = cutset.numCuts();
+  HighsInt origselectednnz = cutset.ARindex_.size();
+  HighsInt selectednnz = origselectednnz;
 
   for (const std::pair<double, HighsInt>& p : efficacious_cuts) {
     bool discard = false;
@@ -415,8 +418,8 @@ void HighsCutPool::separate(const std::vector<double>& sol, HighsDomain& domain,
   assert(int(cutset.ARvalue_.size()) == selectednnz);
   assert(int(cutset.ARindex_.size()) == selectednnz);
 
-  HighsInt offset = 0;
-  for (HighsInt i = 0; i != cutset.numCuts(); ++i) {
+  HighsInt offset = origselectednnz;
+  for (HighsInt i = orignumcuts; i != cutset.numCuts(); ++i) {
     cutset.ARstart_[i] = offset;
     HighsInt cut = cutset.cutindices[i];
     HighsInt start = matrix_.getRowStart(cut);

highs/mip/HighsCutPool.h

@@ -124,7 +124,7 @@ class HighsCutPool {
   double getParallelism(HighsInt row1, HighsInt row2,
                         const HighsCutPool& pool2) const;
 
-  void performAging(bool parallel_sepa = false);
+  void performAging();
 
   void lpCutRemoved(HighsInt cut, bool thread_safe = false);
 

highs/mip/HighsMipSolver.cpp

@@ -143,7 +143,6 @@ void HighsMipSolver::run() {
   HighsMipWorker& master_worker = mipdata_->workers.at(0);
 
 restart:
-  mipdata_->parallel_lock = false;
   if (modelstatus_ == HighsModelStatus::kNotset) {
     // Check limits have not been reached before evaluating root node
     if (mipdata_->checkLimits()) {
@@ -315,7 +314,6 @@ void HighsMipSolver::run() {
   // Initialize worker relaxations and mipworkers
   // todo lps and workers are still empty right now
 
-  mipdata_->parallel_lock = true;
   const HighsInt mip_search_concurrency = options_mip_->mip_search_concurrency;
   const HighsInt num_worker = mip_search_concurrency - 1;
   highs::parallel::TaskGroup tg;
@@ -412,6 +410,14 @@ void HighsMipSolver::run() {
     return search.performed_dive_;
   };
 
+  auto setParallelLock = [&](bool lock) -> void {
+    if (mipdata_->workers.size() <= 1) return;
+    mipdata_->parallel_lock = lock;
+    for (HighsConflictPool& conflictpool : mipdata_->conflictpools) {
+      conflictpool.setAgeLock(lock);
+    }
+  };
+
   auto syncSolutions = [&]() -> void {
     for (HighsMipWorker& worker : mipdata_->workers) {
       for (auto& sol : worker.solutions_) {
@@ -427,6 +433,16 @@ void HighsMipSolver::run() {
     }
   };
 
+  auto syncPools = [&]() -> void {
+    if (mipdata_->workers.size() <= 1 || mipdata_->parallelLockActive()) return;
+    for (HighsInt i = 1; i < mipdata_->conflictpools.size(); ++i) {
+      mipdata_->conflictpools[i].syncConflictPool(mipdata_->conflictPool);
+    }
+    for (HighsInt i = 1; i < mipdata_->cutpools.size(); ++i) {
+      mipdata_->cutpools[i].syncCutPool(*this, mipdata_->cutpool);
+    }
+  };
+
   auto syncGlobalDomain = [&]() -> void {
     if (mipdata_->workers.size() <= 1) return;
     for (HighsMipWorker& worker : mipdata_->workers) {
@@ -515,6 +531,7 @@ void HighsMipSolver::run() {
   auto evaluateNodes = [&](std::vector<HighsInt>& search_indices) -> void {
     std::vector<HighsSearch::NodeResult> search_results(search_indices.size());
     analysis_.mipTimerStart(kMipClockEvaluateNode1);
+    setParallelLock(true);
     for (HighsInt i = 0; i != search_indices.size(); i++) {
       // TODO MT: Remove this dummy if statement
       if (i != 0) continue;
@@ -529,6 +546,7 @@ void HighsMipSolver::run() {
       }
     }
     if (mipdata_->parallelLockActive()) tg.taskWait();
+    setParallelLock(false);
     analysis_.mipTimerStop(kMipClockEvaluateNode1);
     for (HighsInt i = 0; i != search_indices.size(); i++) {
       // TODO MT: Remove this dummy if statement
@@ -629,6 +647,7 @@ void HighsMipSolver::run() {
     std::deque<bool> infeasible;
     std::deque<bool> flush;
     std::vector<bool> prune(search_indices.size(), false);
+    setParallelLock(true);
     for (HighsInt i = 0; i < search_indices.size(); i++) {
       // TODO MT: Remove this redundant if
       if (search_indices[i] != 0 || !mipdata_->workers[search_indices[i]]
@@ -658,6 +677,7 @@ void HighsMipSolver::run() {
         }
       }
     }
+    setParallelLock(false);
 
     // Remove search indices that need a new node
     HighsInt num_search_indices = static_cast<HighsInt>(search_indices.size());
@@ -700,6 +720,7 @@ void HighsMipSolver::run() {
       [&](std::vector<HighsInt>& search_indices) -> bool {
     // the node is still not fathomed, so perform separation
     analysis_.mipTimerStart(kMipClockNodeSearchSeparation);
+    setParallelLock(true);
     for (HighsInt i : search_indices) {
       // TODO MT: Get rid of this line
       if (i != 0) continue;
@@ -715,6 +736,7 @@ void HighsMipSolver::run() {
     }
     analysis_.mipTimerStop(kMipClockNodeSearchSeparation);
     if (mipdata_->parallelLockActive()) tg.taskWait();
+    setParallelLock(false);
 
     for (HighsInt i : search_indices) {
       // TODO MT: Get rid of this line
@@ -765,6 +787,7 @@ void HighsMipSolver::run() {
                                    -analysis_.mipTimerRead(kMipClockTheDive));
     std::vector<HighsSearch::NodeResult> dive_results(
         mip_search_concurrency, HighsSearch::NodeResult::kBranched);
+    setParallelLock(true);
     if (mip_search_concurrency > 1) {
       for (int i = 0; i < mip_search_concurrency; i++) {
         tg.spawn([&, i]() {
@@ -784,6 +807,7 @@ void HighsMipSolver::run() {
         dive_times[0] += analysis_.mipTimerRead(kMipClockNodeSearch);
       }
     }
+    setParallelLock(false);
     bool suboptimal = false;
     for (int i = 0; i < mip_search_concurrency; i++) {
       if (dive_times[i] != -1) {

highs/mip/HighsSeparation.cpp

@@ -140,6 +140,11 @@ HighsInt HighsSeparation::separationRound(HighsDomain& propdomain,
 
   mipworker_.cutpool_->separate(sol.col_value, propdomain, cutset,
                                 mipdata.feastol, mipdata.cutpools);
+  // Also separate the global cut pool
+  if (mipworker_.cutpool_ != &mipdata.cutpool) {
+    mipdata.cutpool.separate(sol.col_value, propdomain, cutset, mipdata.feastol,
+                             mipdata.cutpools, true);
+  }
 
   if (cutset.numCuts() > 0) {
     ncuts += cutset.numCuts();