Forward-merge release/26.02 into main

cpp/include/cuopt/linear_programming/constants.h

@@ -52,12 +52,20 @@
 #define CUOPT_MIP_ABSOLUTE_TOLERANCE          "mip_absolute_tolerance"
 #define CUOPT_MIP_RELATIVE_TOLERANCE          "mip_relative_tolerance"
 #define CUOPT_MIP_INTEGRALITY_TOLERANCE       "mip_integrality_tolerance"
-#define CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING "mip_batch_pdlp_strong_branching"
 #define CUOPT_MIP_ABSOLUTE_GAP                "mip_absolute_gap"
 #define CUOPT_MIP_RELATIVE_GAP                "mip_relative_gap"
 #define CUOPT_MIP_HEURISTICS_ONLY             "mip_heuristics_only"
 #define CUOPT_MIP_SCALING                     "mip_scaling"
 #define CUOPT_MIP_PRESOLVE                    "mip_presolve"
+#define CUOPT_MIP_CUT_PASSES                  "mip_cut_passes"
+#define CUOPT_MIP_MIXED_INTEGER_ROUNDING_CUTS "mip_mixed_integer_rounding_cuts"
+#define CUOPT_MIP_MIXED_INTEGER_GOMORY_CUTS   "mip_mixed_integer_gomory_cuts"
+#define CUOPT_MIP_KNAPSACK_CUTS               "mip_knapsack_cuts"
+#define CUOPT_MIP_STRONG_CHVATAL_GOMORY_CUTS  "mip_strong_chvatal_gomory_cuts"
+#define CUOPT_MIP_REDUCED_COST_STRENGTHENING  "mip_reduced_cost_strengthening"
+#define CUOPT_MIP_CUT_CHANGE_THRESHOLD        "mip_cut_change_threshold"
+#define CUOPT_MIP_CUT_MIN_ORTHOGONALITY       "mip_cut_min_orthogonality"
+#define CUOPT_MIP_BATCH_PDLP_STRONG_BRANCHING "mip_batch_pdlp_strong_branching"
 #define CUOPT_SOLUTION_FILE                   "solution_file"
 #define CUOPT_NUM_CPU_THREADS                 "num_cpu_threads"
 #define CUOPT_NUM_GPUS                        "num_gpus"

cpp/include/cuopt/linear_programming/mip/solver_settings.hpp

@@ -83,12 +83,22 @@ class mip_solver_settings_t {
   friend class problem_checking_t;
   tolerances_t tolerances;
 
-  f_t time_limit       = std::numeric_limits<f_t>::infinity();
-  bool heuristics_only = false;
-  i_t num_cpu_threads  = -1;  // -1 means use default number of threads in branch and bound
-  i_t num_gpus         = 1;
+  f_t time_limit                = std::numeric_limits<f_t>::infinity();
+  i_t node_limit                = std::numeric_limits<i_t>::max();
+  bool heuristics_only          = false;
+  i_t num_cpu_threads           = -1;  // -1 means use default number of threads in branch and bound
+  i_t max_cut_passes            = 10;  // number of cut passes to make
+  i_t mir_cuts                  = -1;
+  i_t mixed_integer_gomory_cuts = -1;
+  i_t knapsack_cuts             = -1;
+  i_t strong_chvatal_gomory_cuts      = -1;
+  i_t reduced_cost_strengthening      = -1;
+  f_t cut_change_threshold            = 1e-3;
+  f_t cut_min_orthogonality           = 0.5;
   i_t mip_batch_pdlp_strong_branching = 0;
+  i_t num_gpus                        = 1;
   bool log_to_console                 = true;
+
   std::string log_file;
   std::string sol_file;
   std::string user_problem_file;

cpp/src/dual_simplex/CMakeLists.txt

@@ -10,6 +10,7 @@ set(DUAL_SIMPLEX_SRC_FILES
   ${CMAKE_CURRENT_SOURCE_DIR}/basis_updates.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/bound_flipping_ratio_test.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/branch_and_bound.cpp
+  ${CMAKE_CURRENT_SOURCE_DIR}/cuts.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/crossover.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/folding.cpp
   ${CMAKE_CURRENT_SOURCE_DIR}/initial_basis.cpp

cpp/src/dual_simplex/barrier.cu

@@ -681,7 +681,7 @@ class iteration_data_t {
           solve_status = chol->solve(U_col, M_col);
           if (solve_status != 0) { return solve_status; }
           if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
-            return -2;
+            return CONCURRENT_HALT_RETURN;
           }
           M.set_column(k, M_col);
 
@@ -700,7 +700,7 @@ class iteration_data_t {
           AD_dense.transpose_multiply(
             1.0, M.values.data() + k * M.m, 0.0, H.values.data() + k * H.m);
           if (settings_.concurrent_halt != nullptr && *settings_.concurrent_halt == 1) {
-            return -2;
+            return CONCURRENT_HALT_RETURN;
           }
         }
 
@@ -1745,7 +1745,7 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   } else {
     status = data.chol->factorize(data.device_ADAT);
   }
-  if (status == -2) { return -2; }
+  if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; }
   if (status != 0) {
     settings.log.printf("Initial factorization failed\n");
     return -1;
@@ -2309,7 +2309,7 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
     data.num_factorizations++;
 
     data.has_solve_info = false;
-    if (status == -2) { return -2; }
+    if (status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; }
 
     if (status < 0) {
       settings.log.printf("Factorization failed.\n");
@@ -2411,7 +2411,7 @@ i_t barrier_solver_t<i_t, f_t>::gpu_compute_search_direction(iteration_data_t<i_
       // TODO Chris, we need to write to cpu because dx is used outside
       // Can't we also GPUify what's usinng this dx?
       raft::copy(dy.data(), data.d_dy_.data(), dy.size(), stream_view_);
-      if (solve_status == -2) { return -2; }
+      if (solve_status == CONCURRENT_HALT_RETURN) { return CONCURRENT_HALT_RETURN; }
       if (solve_status < 0) {
         settings.log.printf("Linear solve failed\n");
         return -1;

cpp/src/dual_simplex/basis_solves.cpp

@@ -363,7 +363,7 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
                                  S_perm_inv);
         if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) {
           settings.log.printf("Concurrent halt\n");
-          return -1;
+          return CONCURRENT_HALT_RETURN;
         }
         if (Srank != Sdim) {
           // Get the rank deficient columns
@@ -582,7 +582,7 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
   }
   if (settings.concurrent_halt != nullptr && *settings.concurrent_halt == 1) {
     settings.log.printf("Concurrent halt\n");
-    return -1;
+    return CONCURRENT_HALT_RETURN;
   }
   if (verbose) {
     printf("Right Lnz+Unz %d t %.3f\n", L.col_start[m] + U.col_start[m], toc(fact_start));
@@ -619,12 +619,13 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
                  const std::vector<i_t>& slacks_needed,
                  std::vector<i_t>& basis_list,
                  std::vector<i_t>& nonbasic_list,
+                 std::vector<i_t>& superbasic_list,
                  std::vector<variable_status_t>& vstatus)
 {
   const i_t m = A.m;
   const i_t n = A.n;
   assert(basis_list.size() == m);
-  assert(nonbasic_list.size() == n - m);
+  assert(nonbasic_list.size() + superbasic_list.size() == n - m);
 
   // Create slack_map
   std::vector<i_t> slack_map(m);  // slack_map[i] = j if column j is e_i
@@ -649,26 +650,39 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
   for (i_t k = 0; k < num_nonbasic; ++k) {
     nonbasic_map[nonbasic_list[k]] = k;
   }
+  // Create a superbasic_map
+  std::vector<i_t> superbasic_map(
+    n, -1);  // superbasic_map[j] = p if superbasic[p] = j, -1 if j is basic/nonbasic
+  const i_t num_superbasic = superbasic_list.size();
+  for (i_t k = 0; k < num_superbasic; ++k) {
+    superbasic_map[superbasic_list[k]] = k;
+  }
 
   const i_t columns_to_replace = deficient.size();
   for (i_t k = 0; k < columns_to_replace; ++k) {
     const i_t bad_j          = basis_list[deficient[k]];
     const i_t replace_i      = slacks_needed[k];
     const i_t replace_j      = slack_map[replace_i];
     basis_list[deficient[k]] = replace_j;
-    assert(nonbasic_map[replace_j] != -1);
-    nonbasic_list[nonbasic_map[replace_j]] = bad_j;
-    vstatus[replace_j]                     = variable_status_t::BASIC;
-    // This is the main issue. What value should bad_j take on.
-    if (lower[bad_j] == -inf && upper[bad_j] == inf) {
-      vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
-    } else if (lower[bad_j] > -inf) {
-      vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
-    } else if (upper[bad_j] < inf) {
-      vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
+    if (nonbasic_map[replace_j] != -1) {
+      nonbasic_list[nonbasic_map[replace_j]] = bad_j;
+      // This is the main issue. What value should bad_j take on.
+      if (lower[bad_j] == -inf && upper[bad_j] == inf) {
+        vstatus[bad_j] = variable_status_t::NONBASIC_FREE;
+      } else if (lower[bad_j] > -inf) {
+        vstatus[bad_j] = variable_status_t::NONBASIC_LOWER;
+      } else if (upper[bad_j] < inf) {
+        vstatus[bad_j] = variable_status_t::NONBASIC_UPPER;
+      } else {
+        assert(1 == 0);
+      }
+    } else if (superbasic_map[replace_j] != -1) {
+      superbasic_list[superbasic_map[replace_j]] = bad_j;
+      vstatus[bad_j]                             = variable_status_t::SUPERBASIC;
     } else {
-      assert(1 == 0);
+      assert(nonbasic_map[replace_j] != -1 || superbasic_map[replace_j] != -1);
     }
+    vstatus[replace_j] = variable_status_t::BASIC;
   }
 
   return 0;
@@ -865,6 +879,7 @@ template int basis_repair<int, double>(const csc_matrix_t<int, double>& A,
                                        const std::vector<int>& slacks_needed,
                                        std::vector<int>& basis_list,
                                        std::vector<int>& nonbasic_list,
+                                       std::vector<int>& superbasic_list,
                                        std::vector<variable_status_t>& vstatus);
 
 template int form_b<int, double>(const csc_matrix_t<int, double>& A,

cpp/src/dual_simplex/basis_solves.hpp

@@ -48,6 +48,7 @@ i_t basis_repair(const csc_matrix_t<i_t, f_t>& A,
                  const std::vector<i_t>& slacks_needed,
                  std::vector<i_t>& basis_list,
                  std::vector<i_t>& nonbasic_list,
+                 std::vector<i_t>& superbasic_list,
                  std::vector<variable_status_t>& vstatus);
 
 // Form the basis matrix B = A(:, basic_list)