diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index 2120152c53..c6f0ec4868 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -565,12 +565,11 @@ i_t factorize_basis(const csc_matrix_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 C(m, m, 1); @@ -591,7 +590,7 @@ i_t factorize_basis(const csc_matrix_t& A, assert(norm_diff < 1e-3); } - return rank; + return (rank == m ? m : -1); } template diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 6e468f9154..c5a5dfe710 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/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(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) { diff --git a/cpp/src/dual_simplex/sparse_matrix.cpp b/cpp/src/dual_simplex/sparse_matrix.cpp index bcd2eb61ea..427f1049d3 100644 --- a/cpp/src/dual_simplex/sparse_matrix.cpp +++ b/cpp/src/dual_simplex/sparse_matrix.cpp @@ -504,18 +504,32 @@ i_t scatter(const csc_matrix_t& A, } template -void csc_matrix_t::check_matrix() const +i_t csc_matrix_t::check_matrix() const { std::vector 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 diff --git a/cpp/src/dual_simplex/sparse_matrix.hpp b/cpp/src/dual_simplex/sparse_matrix.hpp index 6b583e7190..a3db2f8ebe 100644 --- a/cpp/src/dual_simplex/sparse_matrix.hpp +++ b/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;