Fix a bug in basis repair. Recover from numerical issues in primal update. (#249)

This PR includes two fixes:

* Basis repair was failing due to a bug in the right-looking LU where the column permutation was not fully populated when the basis was rank deficient. 

* When computing delta x using a hypersparse solve the solution vector sometimes did not contain a nonzero coefficient associated with the basic leaving index. This was an issue, as we divide by that coefficient to compute the primal step length.  This coefficient may have been missing due to the numerical drop tolerances. To recover from this we perform a regular solve. If the coefficient is small, but nonzero, we can recover. 

In addition this PR adds the option to use threshold rook pivoting inside the right-looking LU factorization. TRP is rank-revealing. However, it produces factors L and U with more fill and is thus about 3% slower. We add TRP in this PR, but keep it disabled behind a flag. A future PR may want to use TRP when basis repair fails. 

New or existing tests cover these changes. The documentation is not affected.

Authors:
  - Chris Maes (https://github.com/chris-maes)

Approvers:
  - Akif ÇÖRDÜK (https://github.com/akifcorduk)

URL: #249

Author: chris-maes

cpp/src/dual_simplex/basis_solves.cpp

@@ -368,6 +368,7 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
           S, settings.threshold_partial_pivoting_tol, identity, S_col_perm, SL, SU, S_perm_inv);
         if (Srank != Sdim) {
           // Get the rank deficient columns
+          deficient.clear();
           deficient.resize(Sdim - Srank);
           for (i_t h = Srank; h < Sdim; ++h) {
             deficient[h - Srank] = col_perm[num_singletons + S_col_perm[h]];

cpp/src/dual_simplex/phase2.cpp

@@ -1604,8 +1604,20 @@ i_t compute_delta_x(const lp_problem_t<i_t, f_t>& lp,
   f_t scale = scaled_delta_xB_sparse.find_coefficient(basic_leaving_index);
   if (scale != scale) {
     // We couldn't find a coefficient for the basic leaving index.
-    // Either this is a bug or the primal step length is inf.
-    return -1;
+    // The coefficient might be very small. Switch to a regular solve and try to recover.
+    std::vector<f_t> rhs;
+    rhs_sparse.to_dense(rhs);
+    const i_t m = basic_list.size();
+    std::vector<f_t> scaled_delta_xB(m);
+    ft.b_solve(rhs, scaled_delta_xB);
+    if (scaled_delta_xB[basic_leaving_index] != 0.0 &&
+        !std::isnan(scaled_delta_xB[basic_leaving_index])) {
+      scaled_delta_xB_sparse.from_dense(scaled_delta_xB);
+      scaled_delta_xB_sparse.negate();
+      scale = -scaled_delta_xB[basic_leaving_index];
+    } else {
+      return -1;
+    }
   }
   const f_t primal_step_length = delta_x_leaving / scale;
   const i_t scaled_delta_xB_nz = scaled_delta_xB_sparse.i.size();
@@ -2856,26 +2868,56 @@ dual::status_t dual_phase2(i_t phase,
     phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 6);
 #endif
     if (should_refactor) {
+      bool should_recompute_x = false;
       if (factorize_basis(lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) ==
           -1) {
+        should_recompute_x = true;
+        settings.log.printf("Failed to factorize basis. Iteration %d\n", iter);
         if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
         basis_repair(lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
         i_t count = 0;
         while (factorize_basis(
                  lp.A, settings, basic_list, L, U, p, pinv, q, deficient, slacks_needed) == -1) {
-          settings.log.printf("Failed to repair basis. %d deficient columns.\n",
+          settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
+                              iter,
                               static_cast<int>(deficient.size()));
           if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
           settings.threshold_partial_pivoting_tol = 1.0;
           count++;
-          if (count > 100) { return dual::status_t::NUMERICAL; }
+          if (count > 10) { return dual::status_t::NUMERICAL; }
           basis_repair(
             lp.A, settings, deficient, slacks_needed, basic_list, nonbasic_list, vstatus);
+
+#ifdef CHECK_BASIS_REPAIR
+          csc_matrix_t<i_t, f_t> B(m, m, 0);
+          form_b(lp.A, basic_list, B);
+          for (i_t k = 0; k < deficient.size(); ++k) {
+            const i_t j         = deficient[k];
+            const i_t col_start = B.col_start[j];
+            const i_t col_end   = B.col_start[j + 1];
+            const i_t col_nz    = col_end - col_start;
+            if (col_nz != 1) {
+              settings.log.printf("Deficient column %d has %d nonzeros\n", j, col_nz);
+            }
+            const i_t i = B.i[col_start];
+            if (i != slacks_needed[k]) {
+              settings.log.printf("Slack %d needed but found %d instead\n", slacks_needed[k], i);
+            }
+          }
+#endif
         }
+
+        settings.log.printf("Successfully repaired basis. Iteration %d\n", iter);
       }
       reorder_basic_list(q, basic_list);
       ft.reset(L, U, p);
       phase2::reset_basis_mark(basic_list, nonbasic_list, basic_mark, nonbasic_mark);
+      if (should_recompute_x) {
+        std::vector<f_t> unperturbed_x(n);
+        phase2::compute_primal_solution_from_basis(
+          lp, ft, basic_list, nonbasic_list, vstatus, unperturbed_x);
+        x = unperturbed_x;
+      }
       phase2::compute_initial_primal_infeasibilities(
         lp, settings, basic_list, x, squared_infeasibilities, infeasibility_indices);
     }