Port of julia #179

src/solver/implementations/default/info.rs

@@ -92,27 +92,29 @@
         let normq = data.get_normq();
 
         // shortcuts for the equilibration matrices
+        let d = &data.equilibration.d;
+        let e = &data.equilibration.e;
         let dinv = &data.equilibration.dinv;
         let einv = &data.equilibration.einv;
-        let cscale = data.equilibration.c;
+        let cinv = T::recip(data.equilibration.c);
 
         // primal and dual costs. dot products are invariant w.r.t
         // equilibration, but we still need to back out the overall
         // objective scaling term c
 
         let xPx_τinvsq_over2 = residuals.dot_xPx * τinv * τinv / (2.).as_T();
-        self.cost_primal = (residuals.dot_qx * τinv + xPx_τinvsq_over2) / cscale;
-        self.cost_dual = (-residuals.dot_bz * τinv - xPx_τinvsq_over2) / cscale;
+        self.cost_primal = (residuals.dot_qx * τinv + xPx_τinvsq_over2) * cinv;
+        self.cost_dual = (-residuals.dot_bz * τinv - xPx_τinvsq_over2) * cinv;
 
         // variables norms, undoing the equilibration.  Do not unscale
         // by τ yet because the infeasibility residuals are ratios of
         // terms that have no affine parts anyway
-        let mut normx = variables.x.norm_scaled(dinv);
-        let mut normz = variables.z.norm_scaled(einv);
+        let mut normx = variables.x.norm_scaled(d);
+        let mut normz = variables.z.norm_scaled(e) * cinv;
         let mut norms = variables.s.norm_scaled(einv);
 
         // primal and dual infeasibility residuals.
-        self.res_primal_inf = residuals.rx_inf.norm_scaled(dinv) / T::max(T::one(), normz);
+        self.res_primal_inf = (residuals.rx_inf.norm_scaled(dinv) * cinv) / T::max(T::one(), normz);
         self.res_dual_inf = T::max(
             residuals.Px.norm_scaled(dinv) / T::max(T::one(), normx),
             residuals.rz_inf.norm_scaled(einv) / T::max(T::one(), normx + norms),
@@ -127,7 +129,7 @@
         self.res_primal =
             residuals.rz.norm_scaled(einv) * τinv / T::max(T::one(), normb + normx + norms);
         self.res_dual =
-            residuals.rx.norm_scaled(dinv) * τinv / T::max(T::one(), normq + normx + normz);
+            residuals.rx.norm_scaled(dinv) * τinv * cinv / T::max(T::one(), normq + normx + normz);
 
         // absolute and relative gaps
         self.gap_abs = T::abs(self.cost_primal - self.cost_dual);
@@ -137,8 +139,8 @@
                 T::min(T::abs(self.cost_primal), T::abs(self.cost_dual)),
             );
 
-        // κ/τ
-        self.ktratio = variables.κ / variables.τ;
+        // κ/τ ratio (scaled)
+        self.ktratio = variables.κ * τinv;
 
         // solve time so far (includes setup)
         self.solve_time = timers.total_time().as_secs_f64();
@@ -169,12 +171,15 @@
             }
 
             // Going backwards. Stop immediately if residuals diverge out of feasibility tolerance.
-            if (self.res_dual > settings.tol_feas
-                && self.res_dual > self.prev_res_dual * (100.).as_T())
-                || (self.res_primal > settings.tol_feas
-                    && self.res_primal > self.prev_res_primal * (100.).as_T())
-            {
-                self.status = SolverStatus::InsufficientProgress;
+            #[allow(clippy::collapsible_if)] // nested if for readability
+            if self.ktratio < T::one() {
+                if (self.res_dual > settings.tol_feas
+                    && self.res_dual > self.prev_res_dual * (100.).as_T())
+                    || (self.res_primal > settings.tol_feas
+                        && self.res_primal > self.prev_res_primal * (100.).as_T())
+                {
+                    self.status = SolverStatus::InsufficientProgress;
+                }
             }
         }
 

src/solver/implementations/default/settings.rs

@@ -53,7 +53,11 @@ pub struct DefaultSettings<T: FloatT> {
     pub tol_ktratio: T,
 
     ///reduced absolute duality gap tolerance
-    #[builder(default = "(5e-5).as_T()")]
+    // NB: reduced_tol_infeas_abs is *smaller* when relaxed, since
+    // we are checking that we are this far into the interior of
+    // an inequality when checking.   Smaller for this value means
+    // "less margin required"
+    #[builder(default = "(5e-12).as_T()")]
     pub reduced_tol_gap_abs: T,
 
     ///reduced relative duality gap tolerance

src/solver/implementations/default/variables.rs

@@ -273,10 +273,10 @@ where
         // also undo the equilibration
         let d = &data.equilibration.d;
         let (e, einv) = (&data.equilibration.e, &data.equilibration.einv);
-        let cscale = data.equilibration.c;
+        let cinv = T::recip(data.equilibration.c);
 
         self.x.hadamard(d).scale(scaleinv);
-        self.z.hadamard(e).scale(scaleinv / cscale);
+        self.z.hadamard(e).scale(scaleinv * cinv);
         self.s.hadamard(einv).scale(scaleinv);
 
         self.τ *= scaleinv;