proxqp: fix absolute value check on duality gap

Author: jcarpent

include/proxsuite/proxqp/dense/solver.hpp

@@ -1073,26 +1073,11 @@ qp_solve( //
       ruiz.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, qpresults.y });
       ruiz.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, qpresults.z });
     }
-    if (is_primal_feasible) {
-
-      if (dual_feasibility_lhs >=
-            qpsettings.refactor_dual_feasibility_threshold &&
-          qpresults.info.rho != qpsettings.refactor_rho_threshold) {
-
-        T rho_new(qpsettings.refactor_rho_threshold);
-
-        refactorize(qpmodel, qpresults, qpwork, rho_new);
-        qpresults.info.rho_updates += 1;
-
-        qpresults.info.rho = rho_new;
-      }
-      if (is_dual_feasible) {
-        if (qpresults.info.duality_gap <=
-            qpsettings.eps_abs +
-              (qpsettings.eps_abs + qpsettings.eps_rel) * rhs_duality_gap) {
-          qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
-          break;
-        }
+    if (is_primal_feasible && is_dual_feasible) {
+      if (std::abs(qpresults.info.duality_gap) <=
+          qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
+        qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
+        break;
       }
     }
     qpresults.info.iter_ext += 1; // We start a new external loop update
@@ -1170,9 +1155,8 @@ qp_solve( //
              std::max(dual_feasibility_rhs_1, qpwork.dual_feasibility_rhs_2)));
 
       if (is_dual_feasible) {
-        if (qpresults.info.duality_gap <=
-            qpsettings.eps_abs +
-              (qpsettings.eps_abs + qpsettings.eps_rel) * rhs_duality_gap) {
+        if (std::abs(qpresults.info.duality_gap) <=
+            qpsettings.eps_abs + qpsettings.eps_rel * rhs_duality_gap) {
           qpresults.info.status = QPSolverOutput::PROXQP_SOLVED;
           break;
         }

include/proxsuite/proxqp/sparse/solver.hpp

@@ -324,11 +324,11 @@ struct PrimalDualGradResult
 /*!
  * Executes the PROXQP algorithm.
  *
- * @param work solver workspace.
- * @param model QP problem model as defined by the user (without any scaling
+ * @param results solver results.
+ * @param data QP problem model as defined by the user (without any scaling
  * performed).
  * @param settings solver settings.
- * @param results solver results.
+ * @param work solver workspace.
  * @param precond preconditioner.
  */
 template<typename T, typename I, typename P>
@@ -781,7 +781,8 @@ qp_solve(Results<T>& results,
       VEG_BIND( // ?
         auto,
         (primal_feasibility_lhs, dual_feasibility_lhs),
-        detail::unscaled_primal_dual_residual(results,
+        detail::unscaled_primal_dual_residual(work,
+                                              results,
                                               primal_residual_eq_scaled,
                                               primal_residual_in_scaled_lo,
                                               primal_residual_in_scaled_up,
@@ -841,9 +842,8 @@ qp_solve(Results<T>& results,
       }
       if (is_primal_feasible(primal_feasibility_lhs) &&
           is_dual_feasible(dual_feasibility_lhs)) {
-        if (results.info.duality_gap <=
-            settings.eps_abs +
-              (settings.eps_rel + settings.eps_abs) * rhs_duality_gap) {
+        if (std::abs(results.info.duality_gap) <=
+            settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
           results.info.pri_res = primal_feasibility_lhs;
           results.info.dua_res = dual_feasibility_lhs;
           results.info.status = QPSolverOutput::PROXQP_SOLVED;
@@ -886,8 +886,10 @@ qp_solve(Results<T>& results,
             LDLT_TEMP_VEC_UNINIT(bool, new_active_constraints, n_in, stack);
             auto rhs = dw;
 
-            work.active_set_low.array() = primal_residual_in_scaled_lo.array() <= 0;
-            work.active_set_up.array() = primal_residual_in_scaled_up.array() >= 0;
+            work.active_set_low.array() =
+              primal_residual_in_scaled_lo.array() <= 0;
+            work.active_set_up.array() =
+              primal_residual_in_scaled_up.array() >= 0;
             new_active_constraints = work.active_set_low || work.active_set_up;
 
             // active set change
@@ -1211,7 +1213,8 @@ qp_solve(Results<T>& results,
       VEG_BIND(
         auto,
         (primal_feasibility_lhs_new, dual_feasibility_lhs_new),
-        detail::unscaled_primal_dual_residual(results,
+        detail::unscaled_primal_dual_residual(work,
+                                              results,
                                               primal_residual_eq_scaled,
                                               primal_residual_in_scaled_lo,
                                               primal_residual_in_scaled_up,
@@ -1232,9 +1235,8 @@ qp_solve(Results<T>& results,
 
       if (is_primal_feasible(primal_feasibility_lhs_new) &&
           is_dual_feasible(dual_feasibility_lhs_new)) {
-        if (results.info.duality_gap <=
-            settings.eps_abs +
-              (settings.eps_abs + settings.eps_rel) * rhs_duality_gap) {
+        if (std::fabs(results.info.duality_gap) <=
+            settings.eps_abs + settings.eps_rel * rhs_duality_gap) {
           results.info.pri_res = primal_feasibility_lhs_new;
           results.info.dua_res = dual_feasibility_lhs_new;
           results.info.status = QPSolverOutput::PROXQP_SOLVED;
@@ -1274,7 +1276,8 @@ qp_solve(Results<T>& results,
       VEG_BIND(
         auto,
         (_, dual_feasibility_lhs_new_2),
-        detail::unscaled_primal_dual_residual(results,
+        detail::unscaled_primal_dual_residual(work,
+                                              results,
                                               primal_residual_eq_scaled,
                                               primal_residual_in_scaled_lo,
                                               primal_residual_in_scaled_up,