Enable parallelism for root node presolve (NVIDIA#371)

Papilo's problem builder was spending a long amount of time building the constraint matrix.
Rewrote that step -> provides ~1.8x speedup on presolve
With TBB there is another ~1.2x additional speedup so ~2x.
[Parallel Presolve.xlsx](https://github.com/user-attachments/files/22145182/Parallel.Presolve.xlsx)

Authors:
  - Hugo Linsenmaier (https://github.com/hlinsen)

Approvers:
  - Rajesh Gandham (https://github.com/rg20)
  - Ramakrishnap (https://github.com/rgsl888prabhu)

URL: NVIDIA#371

Author: hlinsen

ci/build_wheel_libcuopt.sh

@@ -21,8 +21,8 @@ source rapids-init-pip
 package_name="libcuopt"
 package_dir="python/libcuopt"
 
-# Install Boost
-bash ci/utils/install_boost.sh
+# Install Boost and TBB
+bash ci/utils/install_boost_tbb.sh
 
 export SKBUILD_CMAKE_ARGS="-DCUOPT_BUILD_WHEELS=ON;-DDISABLE_DEPRECATION_WARNING=ON"
 

ci/utils/install_boost.sh ci/utils/install_boost_tbb.shci/utils/install_boost.sh renamed to ci/utils/install_boost_tbb.sh

@@ -17,19 +17,19 @@
 
 set -euo pipefail
 
-# Install Boost
+# Install Boost and TBB
 if [ -f /etc/os-release ]; then
     . /etc/os-release
     if [[ "$ID" == "rocky" ]]; then
-        echo "Detected Rocky Linux. Installing Boost via dnf..."
-        dnf install -y boost-devel
+        echo "Detected Rocky Linux. Installing Boost and TBB via dnf..."
+        dnf install -y boost-devel tbb-devel
         if [[ "$(uname -m)" == "x86_64" ]]; then
             dnf install -y gcc-toolset-14-libquadmath-devel
         fi
     elif [[ "$ID" == "ubuntu" ]]; then
-        echo "Detected Ubuntu. Installing Boost via apt..."
+        echo "Detected Ubuntu. Installing Boost and TBB via apt..."
         apt-get update
-        apt-get install -y libboost-dev
+        apt-get install -y libboost-dev libtbb-dev
     else
         echo "Unknown OS: $ID. Please install Boost development libraries manually."
         exit 1

conda/environments/all_cuda-129_arch-aarch64.yaml

@@ -79,6 +79,7 @@ dependencies:
 - sphinxcontrib-openapi
 - sphinxcontrib-websupport
 - sysroot_linux-aarch64==2.28
+- tbb-devel
 - uvicorn==0.34.*
 - zlib
 - pip:

conda/environments/all_cuda-129_arch-x86_64.yaml

@@ -79,6 +79,7 @@ dependencies:
 - sphinxcontrib-openapi
 - sphinxcontrib-websupport
 - sysroot_linux-64==2.28
+- tbb-devel
 - uvicorn==0.34.*
 - zlib
 - pip:

conda/recipes/libcuopt/recipe.yaml

@@ -57,6 +57,7 @@ cache:
       - cuda-version =${{ cuda_version }}
       - cmake ${{ cmake_version }}
       - ninja
+      - tbb-devel
       - zlib
       - bzip2
     host:
@@ -72,6 +73,7 @@ cache:
       - libcusparse-dev
       - cuda-cudart-dev
       - boost
+      - tbb-devel
       - zlib
       - bzip2
 

cpp/CMakeLists.txt

@@ -177,8 +177,8 @@ FetchContent_Declare(
   SYSTEM
 )
 
+find_package(TBB REQUIRED)
 set(BUILD_TESTING OFF CACHE BOOL "Disable test build for papilo")
-set(TBB OFF CACHE BOOL "Disable TBB")
 set(PAPILO_NO_BINARIES ON)
 option(LUSOL "Disable LUSOL" OFF)
 
@@ -254,6 +254,7 @@ target_include_directories(cuopt
 set(CUOPT_PRIVATE_CUDA_LIBS
   CUDA::curand
   CUDA::cusolver
+  TBB::tbb
   OpenMP::OpenMP_CXX)
 
 list(PREPEND CUOPT_PRIVATE_CUDA_LIBS CUDA::cublasLt)

cpp/cmake/thirdparty/FindTBB.cmake

@@ -0,0 +1,90 @@
+# SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION & AFFILIATES. All rights reserved.
+# SPDX-License-Identifier: Apache-2.0
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# FindTBB.cmake - Find TBB (Threading Building Blocks) library
+#
+# This module defines the following variables:
+#   TBB_FOUND        - True if TBB is found
+#   TBB_INCLUDE_DIRS - TBB include directories
+#   TBB_LIBRARIES    - TBB libraries
+#   TBB::tbb         - Imported target for TBB
+
+# Try pkg-config first
+find_package(PkgConfig QUIET)
+if(PkgConfig_FOUND)
+  pkg_check_modules(PC_TBB QUIET tbb)
+endif()
+
+find_path(TBB_INCLUDE_DIR
+  NAMES tbb/tbb.h
+  PATHS
+    ${PC_TBB_INCLUDE_DIRS}
+    /usr/include
+    /usr/local/include
+    /opt/intel/tbb/include
+    /opt/intel/oneapi/tbb/latest/include
+)
+
+find_library(TBB_LIBRARY
+  NAMES tbb
+  PATHS
+    ${PC_TBB_LIBRARY_DIRS}
+    /usr/lib
+    /usr/lib64
+    /usr/local/lib
+    /usr/local/lib64
+    /opt/intel/tbb/lib
+    /opt/intel/oneapi/tbb/latest/lib
+)
+
+find_library(TBB_MALLOC_LIBRARY
+  NAMES tbbmalloc
+  PATHS
+    ${PC_TBB_LIBRARY_DIRS}
+    /usr/lib
+    /usr/lib64
+    /usr/local/lib
+    /usr/local/lib64
+    /opt/intel/tbb/lib
+    /opt/intel/oneapi/tbb/latest/lib
+)
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(TBB
+  REQUIRED_VARS TBB_INCLUDE_DIR TBB_LIBRARY
+)
+
+if(TBB_FOUND AND NOT TARGET TBB::tbb)
+  add_library(TBB::tbb UNKNOWN IMPORTED)
+  set_target_properties(TBB::tbb PROPERTIES
+    IMPORTED_LOCATION "${TBB_LIBRARY}"
+    INTERFACE_INCLUDE_DIRECTORIES "${TBB_INCLUDE_DIR}"
+  )
+
+  if(TBB_MALLOC_LIBRARY)
+    set_target_properties(TBB::tbb PROPERTIES
+      INTERFACE_LINK_LIBRARIES "${TBB_MALLOC_LIBRARY}"
+    )
+  endif()
+
+  # Add compile definitions from pkg-config if available
+  if(PC_TBB_CFLAGS_OTHER)
+    set_target_properties(TBB::tbb PROPERTIES
+      INTERFACE_COMPILE_OPTIONS "${PC_TBB_CFLAGS_OTHER}"
+    )
+  endif()
+endif()
+
+mark_as_advanced(TBB_INCLUDE_DIR TBB_LIBRARY TBB_MALLOC_LIBRARY)

cpp/src/mip/presolve/third_party_presolve.cpp

@@ -19,9 +19,13 @@
 #include <cuopt/logger.hpp>
 #include <mip/mip_constants.hpp>
 #include <mip/presolve/third_party_presolve.hpp>
+#include <utilities/timer.hpp>
 
+#pragma GCC diagnostic push
+#pragma GCC diagnostic ignored "-Wstringop-overflow"  // ignore boost error for pip wheel build
 #include <papilo/core/Presolve.hpp>
 #include <papilo/core/ProblemBuilder.hpp>
+#pragma GCC diagnostic pop
 
 namespace cuopt::linear_programming::detail {
 
@@ -124,18 +128,32 @@ papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>
     builder.setRowRhsAll(h_constr_ub);
   }
 
+  std::vector<papilo::RowFlags> h_row_flags(h_constr_lb.size());
+  std::vector<std::tuple<i_t, i_t, f_t>> h_entries;
   // Add constraints row by row
   for (size_t i = 0; i < h_constr_lb.size(); ++i) {
     // Get row entries
     i_t row_start   = h_offsets[i];
     i_t row_end     = h_offsets[i + 1];
     i_t num_entries = row_end - row_start;
-    builder.addRowEntries(
-      i, num_entries, h_variables.data() + row_start, h_coefficients.data() + row_start);
-    builder.setRowLhsInf(i, h_constr_lb[i] == -std::numeric_limits<f_t>::infinity());
-    builder.setRowRhsInf(i, h_constr_ub[i] == std::numeric_limits<f_t>::infinity());
-    if (h_constr_lb[i] == -std::numeric_limits<f_t>::infinity()) { builder.setRowLhs(i, 0); }
-    if (h_constr_ub[i] == std::numeric_limits<f_t>::infinity()) { builder.setRowRhs(i, 0); }
+    for (size_t j = 0; j < num_entries; ++j) {
+      h_entries.push_back(
+        std::make_tuple(i, h_variables[row_start + j], h_coefficients[row_start + j]));
+    }
+
+    if (h_constr_lb[i] == -std::numeric_limits<f_t>::infinity()) {
+      h_row_flags[i].set(papilo::RowFlag::kLhsInf);
+    } else {
+      h_row_flags[i].unset(papilo::RowFlag::kLhsInf);
+    }
+    if (h_constr_ub[i] == std::numeric_limits<f_t>::infinity()) {
+      h_row_flags[i].set(papilo::RowFlag::kRhsInf);
+    } else {
+      h_row_flags[i].unset(papilo::RowFlag::kRhsInf);
+    }
+
+    if (h_constr_lb[i] == -std::numeric_limits<f_t>::infinity()) { h_constr_lb[i] = 0; }
+    if (h_constr_ub[i] == std::numeric_limits<f_t>::infinity()) { h_constr_ub[i] = 0; }
   }
 
   for (size_t i = 0; i < h_var_lb.size(); ++i) {
@@ -144,7 +162,24 @@ papilo::Problem<f_t> build_papilo_problem(const optimization_problem_t<i_t, f_t>
     if (h_var_lb[i] == -std::numeric_limits<f_t>::infinity()) { builder.setColLb(i, 0); }
     if (h_var_ub[i] == std::numeric_limits<f_t>::infinity()) { builder.setColUb(i, 0); }
   }
-  return builder.build();
+
+  auto problem = builder.build();
+
+  if (h_entries.size()) {
+    auto constexpr const sorted_entries = true;
+    auto csr_storage = papilo::SparseStorage<f_t>(h_entries, num_rows, num_cols, sorted_entries);
+    problem.setConstraintMatrix(csr_storage, h_constr_lb, h_constr_ub, h_row_flags);
+
+    papilo::ConstraintMatrix<f_t>& matrix = problem.getConstraintMatrix();
+    for (int i = 0; i < problem.getNRows(); ++i) {
+      papilo::RowFlags rowFlag = matrix.getRowFlags()[i];
+      if (!rowFlag.test(papilo::RowFlag::kRhsInf) && !rowFlag.test(papilo::RowFlag::kLhsInf) &&
+          matrix.getLeftHandSides()[i] == matrix.getRightHandSides()[i])
+        matrix.getRowFlags()[i].set(papilo::RowFlag::kEquation);
+    }
+  }
+
+  return problem;
 }
 
 template <typename i_t, typename f_t>
@@ -299,14 +334,16 @@ void set_presolve_methods(papilo::Presolve<f_t>& presolver, problem_category_t c
   presolver.addPresolveMethod(uptr(new papilo::Substitution<f_t>()));
 }
 
-template <typename f_t>
+template <typename i_t, typename f_t>
 void set_presolve_options(papilo::Presolve<f_t>& presolver,
                           problem_category_t category,
                           f_t absolute_tolerance,
                           f_t relative_tolerance,
-                          double time_limit)
+                          double time_limit,
+                          i_t num_cpu_threads)
 {
-  presolver.getPresolveOptions().tlim = time_limit;
+  presolver.getPresolveOptions().tlim    = time_limit;
+  presolver.getPresolveOptions().threads = num_cpu_threads;  //  user setting or  0 (automatic)
 }
 
 template <typename i_t, typename f_t>
@@ -315,7 +352,8 @@ std::pair<optimization_problem_t<i_t, f_t>, bool> third_party_presolve_t<i_t, f_
   problem_category_t category,
   f_t absolute_tolerance,
   f_t relative_tolerance,
-  double time_limit)
+  double time_limit,
+  i_t num_cpu_threads)
 {
   cuopt_expects(
     presolve_calls_ == 0, error_type_t::ValidationError, "Presolve can only be called once");
@@ -330,8 +368,8 @@ std::pair<optimization_problem_t<i_t, f_t>, bool> third_party_presolve_t<i_t, f_
 
   papilo::Presolve<f_t> presolver;
   set_presolve_methods<f_t>(presolver, category);
-  set_presolve_options<f_t>(
-    presolver, category, absolute_tolerance, relative_tolerance, time_limit);
+  set_presolve_options<i_t, f_t>(
+    presolver, category, absolute_tolerance, relative_tolerance, time_limit, num_cpu_threads);
 
   // Disable papilo logs
   presolver.setVerbosityLevel(papilo::VerbosityLevel::kQuiet);

cpp/src/mip/presolve/third_party_presolve.hpp

@@ -31,7 +31,8 @@ class third_party_presolve_t {
     problem_category_t category,
     f_t absolute_tolerance,
     f_t relative_tolerance,
-    double time_limit);
+    double time_limit,
+    i_t num_cpu_threads = 0);
 
   void undo(rmm::device_uvector<f_t>& primal_solution,
             rmm::device_uvector<f_t>& dual_solution,

cpp/src/mip/solve.cu

@@ -200,7 +200,8 @@ mip_solution_t<i_t, f_t> solve_mip(optimization_problem_t<i_t, f_t>& op_problem,
                          cuopt::linear_programming::problem_category_t::MIP,
                          settings.tolerances.absolute_tolerance,
                          settings.tolerances.relative_tolerance,
-                         presolve_time_limit);
+                         presolve_time_limit,
+                         settings.num_cpu_threads);
       if (!feasible) {
         return mip_solution_t<i_t, f_t>(mip_termination_status_t::Infeasible,
                                         solver_stats_t<i_t, f_t>{},