Fix an UMA on ns1952667 due to ignoring rank in factorize_basis when no singletons present

cpp/src/dual_simplex/basis_solves.cpp

@@ -565,12 +565,11 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
   }
   q.resize(m);
   f_t fact_start = tic();
-  right_looking_lu(A, medium_tol, basic_list, q, L, U, pinv);
+  rank           = right_looking_lu(A, medium_tol, basic_list, q, L, U, pinv);
   if (verbose) {
     printf("Right Lnz+Unz %d t %.3f\n", L.col_start[m] + U.col_start[m], toc(fact_start));
   }
   inverse_permutation(pinv, p);
-  rank                    = m;
   constexpr bool check_lu = false;
   if (check_lu) {
     csc_matrix_t<i_t, f_t> C(m, m, 1);
@@ -591,7 +590,7 @@ i_t factorize_basis(const csc_matrix_t<i_t, f_t>& A,
     assert(norm_diff < 1e-3);
   }
 
-  return rank;
+  return (rank == m ? m : -1);
 }
 
 template <typename i_t, typename f_t>

cpp/src/dual_simplex/phase2.cpp

@@ -2884,6 +2884,9 @@ dual::status_t dual_phase2(i_t phase,
           settings.log.printf("Failed to repair basis. Iteration %d. %d deficient columns.\n",
                               iter,
                               static_cast<int>(deficient.size()));
+#ifdef CHECK_L_FACTOR
+          if (L.check_matrix() == -1) { settings.log.printf("Bad L after basis repair\n"); }
+#endif
           if (toc(start_time) > settings.time_limit) { return dual::status_t::TIME_LIMIT; }
           settings.threshold_partial_pivoting_tol = 1.0;
           count++;
@@ -2913,6 +2916,9 @@ dual::status_t dual_phase2(i_t phase,
         settings.log.printf("Successfully repaired basis. Iteration %d\n", iter);
       }
       reorder_basic_list(q, basic_list);
+#ifdef CHECK_L_FACTOR
+      if (L.check_matrix() == -1) { settings.log.printf("Bad L factor\n"); }
+#endif
       ft.reset(L, U, p);
       phase2::reset_basis_mark(basic_list, nonbasic_list, basic_mark, nonbasic_mark);
       if (should_recompute_x) {

cpp/src/dual_simplex/sparse_matrix.cpp

@@ -504,18 +504,32 @@ i_t scatter(const csc_matrix_t<i_t, f_t>& A,
 }
 
 template <typename i_t, typename f_t>
-void csc_matrix_t<i_t, f_t>::check_matrix() const
+i_t csc_matrix_t<i_t, f_t>::check_matrix() const
 {
   std::vector<i_t> row_marker(this->m, -1);
   for (i_t j = 0; j < this->n; ++j) {
     const i_t col_start = this->col_start[j];
     const i_t col_end   = this->col_start[j + 1];
+    if (col_start > col_end || col_start > this->col_start[this->n]) {
+      printf("CSC error: column pointers. col start %d col end %d nz %d\n",
+             col_start,
+             col_end,
+             this->col_start[this->n]);
+      return -1;
+    }
     for (i_t p = col_start; p < col_end; ++p) {
       const i_t i = this->i[p];
-      if (row_marker[i] == j) { printf("CSC error: repeated row index %d in column %d\n", i, j); }
+      if (i < 0 || i >= this->m) {
+        printf("CSC error: row index %d not in range [0, %d]\n", i, this->m - 1);
+      }
+      if (row_marker[i] == j) {
+        printf("CSC error: repeated row index %d in column %d\n", i, j);
+        return -1;
+      }
       row_marker[i] = j;
     }
   }
+  return 0;
 }
 
 template <typename i_t, typename f_t>

cpp/src/dual_simplex/sparse_matrix.hpp

@@ -103,7 +103,7 @@ class csc_matrix_t {
   void print_matrix(FILE* fid) const;
 
   // Ensures no repeated row indices within a column
-  void check_matrix() const;
+  i_t check_matrix() const;
 
   // Writes the matrix to a file in Matrix Market format
   void write_matrix_market(FILE* fid) const;