restore residual check in initialization and clean up code related to initial point

Signed-off-by: YUWEN Chen <yuwchen@nvidia.com>

Author: yuwenchen95

cpp/include/cuopt/mathematical_optimization/utilities/internals.hpp

@@ -142,20 +142,5 @@ enum presolver_t : int {
   PSLP    = CUOPT_PRESOLVE_PSLP
 };
 
-/**
- * @brief Barrier primal-dual initial-point strategy.
- *
- * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems.
- * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992).
- * DualLeastSquares: solve augmented or ADAT dual least-squares system.
- * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization).
- */
-enum barrier_dual_initial_point_t : int {
-  Automatic           = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC,
-  LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO,
-  DualLeastSquares    = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES,
-  SedumiMu            = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU
-};
-
 }  // namespace mathematical_optimization
 }  // namespace cuopt

cpp/src/barrier/barrier.cu

@@ -22,7 +22,6 @@
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/solve.hpp>
 
-#include <cuopt/mathematical_optimization/utilities/internals.hpp>
 #include <linear_algebra/sparse_matrix.hpp>
 #include <math_optimization/tic_toc.hpp>
 #include <math_optimization/types.hpp>
@@ -2359,6 +2358,26 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     }
   }
 
+  // Verify A*x = b
+  dense_vector_t<i_t, f_t> init_primal_residual(lp.num_rows);
+  init_primal_residual = lp.rhs;
+  data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual);
+  data.handle_ptr->get_stream().synchronize();
+#ifdef PRINT_INFO
+  settings.log.printf("||b - A * x||: %.16e\n", vector_norm2<i_t, f_t>(init_primal_residual));
+#endif
+
+  if (data.n_upper_bounds > 0) {
+    dense_vector_t<i_t, f_t> init_bound_residual(data.n_upper_bounds);
+    for (i_t k = 0; k < data.n_upper_bounds; k++) {
+      i_t j                  = data.upper_bounds[k];
+      init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j];
+    }
+#ifdef PRINT_INFO
+    settings.log.printf("|| u - w - x||: %e\n", vector_norm2<i_t, f_t>(init_bound_residual));
+#endif
+  }
+
   float64_t epsilon_adjust = 10.0;
   // Push entries into interior of nonnegative orthant and SOC.
   const bool has_soc   = data.has_cones();
@@ -2467,6 +2486,22 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
     data.v.ensure_positive(epsilon_adjust);
   }
 
+  // Verify A'*y + z - E*v  - Q*x = c
+  dense_vector_t<i_t, f_t> init_dual_residual(lp.num_cols);
+  data.z.pairwise_subtract(data.c, init_dual_residual);
+  if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); }
+  data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual);
+  if (data.n_upper_bounds > 0) {
+    for (i_t k = 0; k < data.n_upper_bounds; k++) {
+      i_t j = data.upper_bounds[k];
+      init_dual_residual[j] -= data.v[k];
+    }
+  }
+#ifdef PRINT_INFO
+  settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n",
+                      vector_norm2<i_t, f_t>(init_dual_residual));
+#endif
+
   // Make sure (w, x, v, z) > 0. Skip free variables being handled directly.
   data.w.ensure_positive(epsilon_adjust);
   std::vector<i_t> nonnegative_variables(data.x.size(), 1);
@@ -2487,43 +2522,6 @@ int barrier_solver_t<i_t, f_t>::initial_point(iteration_data_t<i_t, f_t>& data)
   settings.log.printf("min v %e min z %e\n", data.v.minimum(), data.z.minimum());
 #endif
 
-  // Residual checks below reflect the final initial point, after positivity shifts.
-  // Verify A*x = b
-  dense_vector_t<i_t, f_t> init_primal_residual(lp.num_rows);
-  init_primal_residual = lp.rhs;
-  data.cusparse_view_.spmv(1.0, data.x, -1.0, init_primal_residual);
-  data.handle_ptr->get_stream().synchronize();
-#ifdef PRINT_INFO
-  settings.log.printf("||b - A * x||: %.16e\n", vector_norm2<i_t, f_t>(init_primal_residual));
-#endif
-
-  if (data.n_upper_bounds > 0) {
-    dense_vector_t<i_t, f_t> init_bound_residual(data.n_upper_bounds);
-    for (i_t k = 0; k < data.n_upper_bounds; k++) {
-      i_t j                  = data.upper_bounds[k];
-      init_bound_residual[k] = lp.upper[j] - data.w[k] - data.x[j];
-    }
-#ifdef PRINT_INFO
-    settings.log.printf("|| u - w - x||: %e\n", vector_norm2<i_t, f_t>(init_bound_residual));
-#endif
-  }
-
-  // Verify A'*y + z - E*v  - Q*x = c
-  dense_vector_t<i_t, f_t> init_dual_residual(lp.num_cols);
-  data.z.pairwise_subtract(data.c, init_dual_residual);
-  if (data.Q.n > 0) { matrix_vector_multiply(data.Q, -1.0, data.x, 1.0, init_dual_residual); }
-  data.cusparse_view_.transpose_spmv(1.0, data.y, 1.0, init_dual_residual);
-  if (data.n_upper_bounds > 0) {
-    for (i_t k = 0; k < data.n_upper_bounds; k++) {
-      i_t j = data.upper_bounds[k];
-      init_dual_residual[j] -= data.v[k];
-    }
-  }
-#ifdef PRINT_INFO
-  settings.log.printf("||A^T y + z - E*v - Q*x - c ||: %e\n",
-                      vector_norm2<i_t, f_t>(init_dual_residual));
-#endif
-
   return 0;
 }
 

cpp/src/barrier/barrier.hpp

@@ -8,6 +8,7 @@
 
 #include <linear_algebra/dense_vector.hpp>
 
+#include <cuopt/mathematical_optimization/constants.h>
 #include <dual_simplex/presolve.hpp>
 #include <dual_simplex/simplex_solver_settings.hpp>
 #include <dual_simplex/solution.hpp>
@@ -18,6 +19,21 @@
 #include <rmm/device_uvector.hpp>
 namespace cuopt::mathematical_optimization::barrier {
 
+/**
+ * @brief Barrier primal-dual initial-point strategy.
+ *
+ * Automatic: use Lustig-Marsten-Shanno for LP/QP; Sturm/SeDuMi mu-based point for conic problems.
+ * LustigMarstenShanno: Mehrotra-style dual start (Lustig, Marsten, Shanno, SIAM J. Optim. 1992).
+ * DualLeastSquares: solve augmented or ADAT dual least-squares system.
+ * SedumiMu: Sturm/SeDuMi mu-based primal+dual point (no factorization).
+ */
+enum barrier_dual_initial_point_t : int {
+  Automatic           = CUOPT_BARRIER_DUAL_INITIAL_POINT_AUTOMATIC,
+  LustigMarstenShanno = CUOPT_BARRIER_DUAL_INITIAL_POINT_LUSTIG_MARSTEN_SHANNO,
+  DualLeastSquares    = CUOPT_BARRIER_DUAL_INITIAL_POINT_DUAL_LEAST_SQUARES,
+  SedumiMu            = CUOPT_BARRIER_DUAL_INITIAL_POINT_SEDUMI_MU
+};
+
 /** Validates SOC layout on an simplex::lp_problem_t before barrier presolve/solve. */
 template <typename i_t, typename f_t>
 bool validate_barrier_cone_layout(const simplex::lp_problem_t<i_t, f_t>& problem,

cpp/src/dual_simplex/simplex_solver_settings.hpp

@@ -167,7 +167,8 @@ struct simplex_solver_settings_t {
   i_t augmented;  // -1 automatic, 0 to solve with ADAT, 1 to solve with augmented system
   i_t dualize;    // -1 automatic, 0 to not dualize, 1 to dualize
   i_t ordering;   // -1 automatic, 0 to use nested dissection, 1 to use AMD
-  i_t barrier_dual_initial_point;  // barrier_dual_initial_point_t; see internals.hpp
+  i_t barrier_dual_initial_point;  // -1 automatic, 0 Lustig-Marsten-Shanno,
+                                   // 1 dual least squares, 2 SeDuMi mu-based
   bool check_Q;                    // true to check if Q is positive semidefinite
   bool crossover;                  // true to do crossover, false to not
   i_t refactor_frequency;          // number of basis updates before refactorization