diff --git a/cpp/src/dual_simplex/basis_solves.cpp b/cpp/src/dual_simplex/basis_solves.cpp index 53464fc66..a2243f540 100644 --- a/cpp/src/dual_simplex/basis_solves.cpp +++ b/cpp/src/dual_simplex/basis_solves.cpp @@ -368,6 +368,7 @@ i_t factorize_basis(const csc_matrix_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]]; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index a83e3bf72..76f4768ab 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -1604,8 +1604,20 @@ i_t compute_delta_x(const lp_problem_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 rhs; + rhs_sparse.to_dense(rhs); + const i_t m = basic_list.size(); + std::vector 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(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 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 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); } diff --git a/cpp/src/dual_simplex/right_looking_lu.cpp b/cpp/src/dual_simplex/right_looking_lu.cpp index 7fe276ec6..caf80ad11 100644 --- a/cpp/src/dual_simplex/right_looking_lu.cpp +++ b/cpp/src/dual_simplex/right_looking_lu.cpp @@ -31,12 +31,11 @@ namespace { // submatrix during the LU factorization template struct element_t { - i_t i; // row index - i_t j; // column index - f_t x; // coefficient value - i_t - next_in_column; // index of the next element in the column: nullptr if there is no next element - i_t next_in_row; // index of the next element in the row: nullptr if there is no next element + i_t i; // row index + i_t j; // column index + f_t x; // coefficient value + i_t next_in_column; // index of the next element in the column: kNone if there is no next element + i_t next_in_row; // index of the next element in the row: kNone if there is no next element }; constexpr int kNone = -1; @@ -165,6 +164,34 @@ void initialize_max_in_column(const std::vector& first_in_col, } } +template +f_t maximum_in_row(i_t i, + const std::vector& first_in_row, + std::vector>& elements) +{ + f_t max_in_row = 0.0; + for (i_t p = first_in_row[i]; p != kNone; p = elements[p].next_in_row) { + element_t* entry = &elements[p]; + assert(entry->i == i); + max_in_row = std::max(max_in_row, std::abs(entry->x)); + } + return max_in_row; +} + +template +void initialize_max_in_row(const std::vector& first_in_row, + std::vector>& elements, + std::vector& max_in_row) +{ + const i_t m = first_in_row.size(); + for (i_t i = 0; i < m; ++i) { + max_in_row[i] = maximum_in_row(i, first_in_row, elements); + } +} + +#undef THRESHOLD_ROOK_PIVOTING // Disable threshold rook pivoting for now. + // 3% slower when enabled. But keep it around + // for challenging numerical problems. template i_t markowitz_search(const std::vector& Cdegree, const std::vector& Rdegree, @@ -173,6 +200,7 @@ i_t markowitz_search(const std::vector& Cdegree, const std::vector& first_in_row, const std::vector& first_in_col, const std::vector& max_in_column, + const std::vector& max_in_row, std::vector>& elements, f_t pivot_tol, f_t threshold_tol, @@ -199,6 +227,7 @@ i_t markowitz_search(const std::vector& Cdegree, element_t* entry = &elements[p]; const i_t i = entry->i; assert(entry->j == j); +#ifdef CHECK_RDEGREE if (Rdegree[i] < 0) { if (verbose) { printf("Rdegree[%d] %d. Searching in column %d. Entry i %d j %d val %e\n", @@ -210,9 +239,13 @@ i_t markowitz_search(const std::vector& Cdegree, entry->x); } } +#endif assert(Rdegree[i] >= 0); const i_t Mij = (Rdegree[i] - 1) * (nz - 1); if (Mij < markowitz && std::abs(entry->x) >= threshold_tol * max_in_col && +#ifdef THRESHOLD_ROOK_PIVOTING + std::abs(entry->x) >= threshold_tol * max_in_row[i] && +#endif std::abs(entry->x) >= pivot_tol) { markowitz = Mij; pivot_i = i; @@ -233,6 +266,9 @@ i_t markowitz_search(const std::vector& Cdegree, assert(row_count[nz].size() >= 0); for (const i_t i : row_count[nz]) { assert(Rdegree[i] == nz); +#ifdef THRESHOLD_ROOK_PIVOTING + const f_t max_in_row_i = max_in_row[i]; +#endif for (i_t p = first_in_row[i]; p != kNone; p = elements[p].next_in_row) { element_t* entry = &elements[p]; const i_t j = entry->j; @@ -241,6 +277,9 @@ i_t markowitz_search(const std::vector& Cdegree, assert(Cdegree[j] >= 0); const i_t Mij = (nz - 1) * (Cdegree[j] - 1); if (Mij < markowitz && std::abs(entry->x) >= threshold_tol * max_in_col && +#ifdef THRESHOLD_ROOK_PIVOTING + std::abs(entry->x) >= threshold_tol * max_in_row_i && +#endif std::abs(entry->x) >= pivot_tol) { markowitz = Mij; pivot_i = i; @@ -257,7 +296,7 @@ i_t markowitz_search(const std::vector& Cdegree, nz++; } if (nsearch > 10) { - if (verbose) { printf("nsearch %d\n", nsearch); } + if constexpr (verbose) { printf("nsearch %d\n", nsearch); } } return nsearch; } @@ -333,6 +372,7 @@ void schur_complement(i_t pivot_i, std::vector& row_last_workspace, std::vector& column_j_workspace, std::vector& max_in_column, + std::vector& max_in_row, std::vector& Rdegree, std::vector& Cdegree, std::vector>& row_count, @@ -378,6 +418,9 @@ void schur_complement(i_t pivot_i, e2->x -= val; const f_t abs_e2x = std::abs(e2->x); if (abs_e2x > max_in_column[j]) { max_in_column[j] = abs_e2x; } +#ifdef THRESHOLD_ROOK_PIVOTING + if (abs_e2x > max_in_row[i]) { max_in_row[i] = abs_e2x; } +#endif } else { element_t fill; fill.i = i; @@ -385,6 +428,9 @@ void schur_complement(i_t pivot_i, fill.x = -val; const f_t abs_fillx = std::abs(fill.x); if (abs_fillx > max_in_column[j]) { max_in_column[j] = abs_fillx; } +#ifdef THRESHOLD_ROOK_PIVOTING + if (abs_fillx > max_in_row[i]) { max_in_row[i] = abs_fillx; } +#endif fill.next_in_column = kNone; fill.next_in_row = kNone; elements.push_back(fill); @@ -484,7 +530,7 @@ void remove_pivot_col(i_t pivot_i, i_t pivot_j, std::vector& first_in_col, std::vector& first_in_row, - std::vector& max_in_column, + std::vector& max_in_row, std::vector>& elements) { // Remove the pivot col @@ -492,6 +538,9 @@ void remove_pivot_col(i_t pivot_i, element_t* e = &elements[p1]; const i_t i = e->i; i_t last = kNone; +#ifdef THRESHOLD_ROOK_PIVOTING + f_t max_in_row_i = 0.0; +#endif for (i_t p = first_in_row[i]; p != kNone; p = elements[p].next_in_row) { element_t* entry = &elements[p]; if (entry->j == pivot_j) { @@ -504,8 +553,17 @@ void remove_pivot_col(i_t pivot_i, entry->j = -1; entry->x = std::numeric_limits::quiet_NaN(); } +#ifdef THRESHOLD_ROOK_PIVOTING + else { + const f_t abs_entryx = std::abs(entry->x); + if (abs_entryx > max_in_row_i) { max_in_row_i = abs_entryx; } + } +#endif last = p; } +#ifdef THRESHOLD_ROOK_PIVOTING + max_in_row[i] = max_in_row_i; +#endif } first_in_col[pivot_j] = kNone; } @@ -549,7 +607,11 @@ i_t right_looking_lu(const csc_matrix_t& A, std::vector column_j_workspace(n, kNone); std::vector row_last_workspace(n); std::vector max_in_column(n); + std::vector max_in_row(m); initialize_max_in_column(first_in_col, elements, max_in_column); +#ifdef THRESHOLD_ROOK_PIVOTING + initialize_max_in_row(first_in_row, elements, max_in_row); +#endif csr_matrix_t Urow; // We will store U by rows in Urow during the factorization and // translate back to U at the end @@ -561,10 +623,10 @@ i_t right_looking_lu(const csc_matrix_t& A, L.x.clear(); L.i.clear(); - for (i_t k = 0; k < n; ++k) { - pinv[k] = -1; - q[k] = -1; - } + std::fill(q.begin(), q.end(), -1); + std::fill(pinv.begin(), pinv.end(), -1); + std::vector qinv(n); + std::fill(qinv.begin(), qinv.end(), -1); i_t pivots = 0; for (i_t k = 0; k < n; ++k) { @@ -584,6 +646,7 @@ i_t right_looking_lu(const csc_matrix_t& A, first_in_row, first_in_col, max_in_column, + max_in_row, elements, pivot_tol, threshold_tol, @@ -598,6 +661,7 @@ i_t right_looking_lu(const csc_matrix_t& A, // Pivot pinv[pivot_i] = k; // pivot_i is the kth pivot row q[k] = pivot_j; + qinv[pivot_j] = k; const f_t pivot_val = pivot_entry->x; assert(std::abs(pivot_val) >= pivot_tol); pivots++; @@ -656,6 +720,7 @@ i_t right_looking_lu(const csc_matrix_t& A, row_last_workspace, column_j_workspace, max_in_column, + max_in_row, Rdegree, Cdegree, row_count, @@ -664,7 +729,7 @@ i_t right_looking_lu(const csc_matrix_t& A, // Remove the pivot row remove_pivot_row(pivot_i, pivot_j, first_in_col, first_in_row, max_in_column, elements); - remove_pivot_col(pivot_i, pivot_j, first_in_col, first_in_row, max_in_column, elements); + remove_pivot_col(pivot_i, pivot_j, first_in_col, first_in_row, max_in_row, elements); // Set pivot entry to sentinel value pivot_entry->i = -1; @@ -695,6 +760,30 @@ i_t right_looking_lu(const csc_matrix_t& A, } #endif +#ifdef CHECK_MAX_IN_ROW + // Check that maximum in row is maintained + for (i_t i = 0; i < m; ++i) { + if (Rdegree[i] == -1) { continue; } + const f_t max_in_row_i = max_in_row[i]; + bool found_max = false; + f_t largest_abs_x = 0.0; + for (i_t p = first_in_row[i]; p != kNone; p = elements[p].next_in_row) { + const f_t abs_e2x = std::abs(elements[p].x); + if (abs_e2x > largest_abs_x) { largest_abs_x = abs_e2x; } + if (abs_e2x > max_in_row_i) { + printf("Found max in row %d is %e but %e\n", i, max_in_row_i, abs_e2x); + } + assert(abs_e2x <= max_in_row_i); + if (abs_e2x == max_in_row_i) { found_max = true; } + } + if (!found_max) { + printf( + "Did not find max %e in row %d. Largest abs x is %e\n", max_in_row_i, i, largest_abs_x); + } + assert(found_max); + } +#endif + #if CHECK_BAD_ENTRIES for (Int j = 0; j < n; j++) { for (Int p = first_in_col[j]; p != kNone; p = elements[p].next_in_column) { @@ -761,6 +850,15 @@ i_t right_looking_lu(const csc_matrix_t& A, for (i_t i = 0; i < m; ++i) { if (pinv[i] == -1) { pinv[i] = start++; } } + + // Finalize the permutation q. Do this by first completing the inverse permutation qinv. + // Then invert qinv to get the final permutation q. + start = pivots; + for (i_t j = 0; j < n; ++j) { + if (qinv[j] == -1) { qinv[j] = start++; } + } + inverse_permutation(qinv, q); + return pivots; } @@ -852,7 +950,11 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, std::vector column_j_workspace(m, kNone); std::vector row_last_workspace(m); std::vector max_in_column(n); + std::vector max_in_row(m); initialize_max_in_column(first_in_col, elements, max_in_column); +#ifdef THRESHOLD_ROOK_PIVOTING + initialize_max_in_row(first_in_row, elements, max_in_row); +#endif settings.log.debug("Empty rows %ld\n", row_count[0].size()); settings.log.debug("Empty cols %ld\n", col_count[0].size()); @@ -884,6 +986,7 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, first_in_row, first_in_col, max_in_column, + max_in_row, elements, pivot_tol, threshold_tol, @@ -924,6 +1027,7 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, row_last_workspace, column_j_workspace, max_in_column, + max_in_row, Rdegree, Cdegree, row_count, @@ -933,8 +1037,7 @@ i_t right_looking_lu_row_permutation_only(const csc_matrix_t& A, // Remove the pivot row remove_pivot_row( pivot_i, pivot_j, first_in_col, first_in_row, max_in_column, elements); - remove_pivot_col( - pivot_i, pivot_j, first_in_col, first_in_row, max_in_column, elements); + remove_pivot_col(pivot_i, pivot_j, first_in_col, first_in_row, max_in_row, elements); // Set pivot entry to sentinel value pivot_entry->i = -1;