Crossover dualize fix

cpp/src/dual_simplex/presolve.cpp

@@ -503,8 +503,6 @@ i_t convert_less_than_to_equal(const user_problem_t<i_t, f_t>& user_problem,
   // We must convert rows in the form: a_i^T x <= beta
   // into: a_i^T x + s_i = beta, s_i >= 0
 
-  csr_matrix_t<i_t, f_t> Arow(0, 0, 0);
-  problem.A.to_compressed_row(Arow);
   i_t num_cols = problem.num_cols + less_rows;
   i_t nnz      = problem.A.col_start[problem.num_cols] + less_rows;
   problem.A.col_start.resize(num_cols + 1);

cpp/src/dual_simplex/solve.cpp

@@ -388,18 +388,62 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
       settings.log.printf("Primal objective: %e\n",
                           dot<i_t, f_t>(dualize_info.primal_problem.objective, primal_solution.x));
 
-      std::vector<f_t> primal_residual = dualize_info.primal_problem.rhs;
-      matrix_vector_multiply(
-        dualize_info.primal_problem.A, 1.0, primal_solution.x, -1.0, primal_residual);
       std::vector<i_t> inequality_rows(dualize_info.primal_problem.num_rows, 1);
       for (i_t i : dualize_info.equality_rows) {
         inequality_rows[i] = 0;
       }
+      i_t less_rows = 0;
       for (i_t i = 0; i < dualize_info.primal_problem.num_rows; ++i) {
         if (inequality_rows[i] == 1) {
-          primal_residual[i] = std::max(primal_residual[i], 0.0);  // a_i^T x - b_i <= 0
+          less_rows++;
         }
       }
+      // Add slack variables to the primal problem
+      if (less_rows > 0)
+      {
+        std::vector<f_t> slack_info = dualize_info.primal_problem.rhs;
+        matrix_vector_multiply(dualize_info.primal_problem.A, -1.0, primal_solution.x, 1.0, slack_info);
+
+
+        lp_problem_t<i_t, f_t>& problem = dualize_info.primal_problem;
+        i_t num_cols = problem.num_cols + less_rows;
+        i_t nnz      = problem.A.col_start[problem.num_cols] + less_rows;
+        problem.A.col_start.resize(num_cols + 1);
+        problem.A.i.resize(nnz);
+        problem.A.x.resize(nnz);
+        problem.lower.resize(num_cols);
+        problem.upper.resize(num_cols);
+        problem.objective.resize(num_cols);
+        primal_solution.x.resize(num_cols);
+        primal_solution.z.resize(num_cols);
+
+        i_t p = problem.A.col_start[problem.num_cols];
+        i_t j = problem.num_cols;
+        for (i_t i = 0; i < problem.num_rows; i++) {
+          if (inequality_rows[i] == 1) {
+            problem.lower[j]     = 0.0;
+            problem.upper[j]     = INFINITY;
+            problem.objective[j] = 0.0;
+            problem.A.i[p]       = i;
+            problem.A.x[p]       = 1.0;
+            primal_solution.x[j]  = slack_info[i];
+            primal_solution.z[j]  = -primal_solution.y[i];
+            problem.A.col_start[j++] = p++;
+            inequality_rows[i]   = 0;
+            less_rows--;
+          }
+        }
+        problem.A.col_start[num_cols] = p;
+        assert(less_rows == 0);
+        assert(p == nnz);
+        problem.A.n      = num_cols;
+        problem.num_cols = num_cols;
+      }
+
+      std::vector<f_t> primal_residual = dualize_info.primal_problem.rhs;
+      matrix_vector_multiply(
+        dualize_info.primal_problem.A, 1.0, primal_solution.x, -1.0, primal_residual);
+
       f_t primal_residual_norm     = vector_norm_inf<i_t, f_t>(primal_residual);
       const f_t norm_b             = vector_norm_inf<i_t, f_t>(dualize_info.primal_problem.rhs);
       f_t primal_relative_residual = primal_residual_norm / (1.0 + norm_b);
@@ -436,6 +480,13 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
   if (!settings.crossover) { return barrier_status; }
 
   if (settings.crossover && barrier_status == lp_status_t::OPTIMAL) {
+
+    {
+      std::vector<f_t> rhs = original_lp.rhs;
+      matrix_vector_multiply(original_lp.A, 1.0, lp_solution.x, -1.0, rhs);
+      f_t primal_residual = vector_norm_inf<i_t, f_t>(rhs);
+      settings.log.printf("Primal residual before adding artificial variables: %e\n", primal_residual);
+    }
     // Check to see if we need to add artifical variables
     std::vector<i_t> artificial_variables;
     artificial_variables.reserve(original_lp.num_rows);
@@ -481,6 +532,12 @@ lp_status_t solve_linear_program_with_barrier(const user_problem_t<i_t, f_t>& us
              lp_solution.x.size(),
              lp_solution.z.size());
 #endif
+
+        std::vector<f_t> rhs = original_lp.rhs;
+        matrix_vector_multiply(original_lp.A, 1.0, lp_solution.x, -1.0, rhs);
+        f_t primal_residual = vector_norm_inf<i_t, f_t>(rhs);
+        settings.log.printf("Primal residual after adding artificial variables: %e\n", primal_residual);
+
     }
 
     // Run crossover

cpp/src/mip/diversity/weights.cuh

@@ -17,6 +17,7 @@
 
 #pragma once
 
+#include <thrust/fill.h>
 #include <raft/core/handle.hpp>
 #include <rmm/device_scalar.hpp>
 #include <rmm/device_uvector.hpp>

cpp/src/mip/problem/presolve_data.cuh

@@ -20,6 +20,7 @@
 #include <cuopt/linear_programming/optimization_problem.hpp>
 
 #include <thrust/sequence.h>
+#include <thrust/uninitialized_fill.h>
 #include <rmm/device_uvector.hpp>
 
 namespace cuopt {

cpp/src/mip/relaxed_lp/lp_state.cuh

@@ -17,6 +17,7 @@
 
 #pragma once
 
+#include <thrust/fill.h>
 #include <raft/util/cudart_utils.hpp>
 #include <rmm/device_uvector.hpp>