diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 33a2d983c..dee39ca50 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -2506,7 +2506,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut set_uninitialized_steepest_edge_norms(original_lp_, basic_list, edge_norms_); pc_.resize(original_lp_.num_cols); - original_lp_.A.transpose(pc_.AT); + pc_.Arow = Arow_; { raft::common::nvtx::range scope_sb("BB::strong_branching"); strong_branching(original_lp_, diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index c38e98e27..aa059b532 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -66,14 +66,14 @@ template objective_change_estimate_t single_pivot_objective_change_estimate( const lp_problem_t& lp, const simplex_solver_settings_t& settings, - const csc_matrix_t& A_transpose, + const csr_matrix_t& Arow, const std::vector& vstatus, i_t variable_j, i_t basic_j, const lp_solution_t& lp_solution, const std::vector& basic_list, const std::vector& nonbasic_list, - const std::vector& nonbasic_mark, + const std::vector& nonbasic_end, basis_update_mpf_t& basis_factors, std::vector& workspace, std::vector& delta_z, @@ -82,7 +82,6 @@ objective_change_estimate_t single_pivot_objective_change_estimate( // Compute the objective estimate for the down and up branches of variable j assert(variable_j >= 0); assert(basic_j >= 0); - // Down branch i_t direction = -1; sparse_vector_t e_k(lp.num_rows, 0); @@ -104,11 +103,11 @@ objective_change_estimate_t single_pivot_objective_change_estimate( std::vector delta_z_indices; // delta_z starts out all zero if (use_transpose) { - compute_delta_z(A_transpose, + compute_delta_z(Arow, delta_y, variable_j, direction, - nonbasic_mark, + nonbasic_end, workspace, delta_z_indices, delta_z, @@ -234,10 +233,8 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, basic_map[basic_list[i]] = i; } - std::vector nonbasic_mark(n, -1); - for (i_t i = 0; i < n - m; i++) { - nonbasic_mark[nonbasic_list[i]] = i; - } + std::vector nonbasic_end(lp.num_rows); + compute_initial_nonbasic_end(basic_map, pc.Arow, nonbasic_end); for (i_t k = 0; k < fractional.size(); k++) { const i_t j = fractional[k]; @@ -246,14 +243,14 @@ void initialize_pseudo_costs_with_estimate(const lp_problem_t& lp, objective_change_estimate_t estimate = single_pivot_objective_change_estimate(lp, settings, - pc.AT, + pc.Arow, vstatus, j, basic_map[j], lp_solution, basic_list, nonbasic_list, - nonbasic_mark, + nonbasic_end, basis_factors, workspace, delta_z, @@ -1476,10 +1473,13 @@ i_t pseudo_costs_t::reliable_variable_selection( basic_map[worker->basic_list[i]] = i; } - std::vector nonbasic_mark(n, -1); - for (i_t i = 0; i < n - m; i++) { - nonbasic_mark[worker->nonbasic_list[i]] = i; - } + // Each thread will have a different basis + // So we need to make a copy of Arow before we permute the columns + // so that nonbasic variables appear first + csr_matrix_t local_Arow = Arow; + + std::vector nonbasic_end(m); + compute_initial_nonbasic_end(basic_map, local_Arow, nonbasic_end); for (auto& [score, j] : unreliable_list) { if (pseudo_cost_num_down[j] == 0 || pseudo_cost_num_up[j] == 0) { @@ -1487,14 +1487,14 @@ i_t pseudo_costs_t::reliable_variable_selection( objective_change_estimate_t estimate = single_pivot_objective_change_estimate(worker->leaf_problem, settings, - AT, + local_Arow, node_ptr->vstatus, j, basic_map[j], leaf_solution, worker->basic_list, worker->nonbasic_list, - nonbasic_mark, + nonbasic_end, worker->basis_factors, workspace, delta_z, diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 009bd8b81..43d6448e6 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -435,7 +435,7 @@ class pseudo_costs_t { pseudo_cost_num_up(num_variables), pseudo_cost_mutex_up(num_variables), pseudo_cost_mutex_down(num_variables), - AT(1, 1, 1) + Arow(1, 1, 1) { } @@ -526,7 +526,7 @@ class pseudo_costs_t { reliability_branching_settings_t reliability_branching_settings; - csc_matrix_t AT; // Transpose of the constraint matrix A + csr_matrix_t Arow; std::vector> pseudo_cost_sum_up; std::vector> pseudo_cost_sum_down; std::vector> pseudo_cost_num_up; diff --git a/cpp/src/dual_simplex/phase2.cpp b/cpp/src/dual_simplex/phase2.cpp index 5b1130796..43be0b881 100644 --- a/cpp/src/dual_simplex/phase2.cpp +++ b/cpp/src/dual_simplex/phase2.cpp @@ -130,79 +130,136 @@ void compute_reduced_cost_update(const lp_problem_t& lp, } template -void compute_delta_z(const csc_matrix_t& A_transpose, +void compute_delta_z(const csr_matrix_t& Arow, const sparse_vector_t& delta_y, i_t leaving_index, i_t direction, - const std::vector& nonbasic_mark, + const std::vector& nonbasic_end, std::vector& delta_z_mark, std::vector& delta_z_indices, std::vector& delta_z, f_t& work_estimate) { // delta_zN = - N'*delta_y - const i_t nz_delta_y = delta_y.i.size(); - size_t nnz_processed = 0; - size_t nonbasic_marked = 0; + const i_t nz_delta_y = delta_y.i.size(); + size_t nnz_processed = 0; for (i_t k = 0; k < nz_delta_y; k++) { const i_t i = delta_y.i[k]; const f_t delta_y_i = delta_y.x[k]; if (std::abs(delta_y_i) < 1e-12) { continue; } - const i_t row_start = A_transpose.col_start[i]; - const i_t row_end = A_transpose.col_start[i + 1]; + const i_t row_start = Arow.row_start[i]; + const i_t row_end = nonbasic_end[i] + 1; nnz_processed += row_end - row_start; for (i_t p = row_start; p < row_end; ++p) { - const i_t j = A_transpose.i[p]; - if (nonbasic_mark[j] >= 0) { - delta_z[j] -= delta_y_i * A_transpose.x[p]; - nonbasic_marked++; - if (!delta_z_mark[j]) { - delta_z_mark[j] = 1; - delta_z_indices.push_back(j); - } + const i_t j = Arow.j[p]; + delta_z[j] -= delta_y_i * Arow.x[p]; + if (!delta_z_mark[j]) { + delta_z_mark[j] = 1; + delta_z_indices.push_back(j); } } } work_estimate += 4 * nz_delta_y; - work_estimate += 2 * nnz_processed; - work_estimate += 3 * nonbasic_marked; + work_estimate += 4 * nnz_processed; work_estimate += 2 * delta_z_indices.size(); // delta_zB = sigma*ei delta_z[leaving_index] = direction; +} -#ifdef CHECK_CHANGE_IN_REDUCED_COST - const i_t m = A_transpose.n; - const i_t n = A_transpose.m; - std::vector delta_y_dense(m); - delta_y.to_dense(delta_y_dense); - std::vector delta_z_check(n); - std::vector delta_z_mark_check(n, 0); - std::vector delta_z_indices_check; - phase2::compute_reduced_cost_update(lp, - basic_list, - nonbasic_list, - delta_y_dense, - leaving_index, - direction, - delta_z_mark_check, - delta_z_indices_check, - delta_z_check, - work_estimate); - f_t error_check = 0.0; - for (i_t k = 0; k < n; ++k) { - const f_t diff = std::abs(delta_z[k] - delta_z_check[k]); - if (diff > 1e-6) { - printf("delta_z error %d transpose %e no transpose %e diff %e\n", - k, - delta_z[k], - delta_z_check[k], - diff); +template +void compute_initial_nonbasic_end(const std::vector& basic_mark, + csr_matrix_t& Arow, + std::vector& nonbasic_end) +{ + i_t m = Arow.m; + // Swap coefficients so that all coefficients for nonbasic variables in row i + // appear from Arow.row_start[i] to nonbasic_end[i] + for (i_t i = 0; i < m; i++) { + i_t lo = Arow.row_start[i]; + i_t hi = Arow.row_start[i + 1] - 1; + while (lo <= hi) { + const i_t j = Arow.j[lo]; + if (basic_mark[j] >= 0) { + // Swap coefficients + std::swap(Arow.j[lo], Arow.j[hi]); + std::swap(Arow.x[lo], Arow.x[hi]); + hi--; + } else { + lo++; + } } - error_check = std::max(error_check, diff); + nonbasic_end[i] = hi; } - if (error_check > 1e-6) { printf("delta_z error %e\n", error_check); } -#endif +} + +template +void update_Arow(i_t leaving, + i_t entering, + const csc_matrix_t& A, + std::vector& row_mark, + std::vector& nonbasic_end, + csr_matrix_t& Arow, + f_t& work_estimate) +{ + // Update the Arow matrix to reflect the new basis. So + // that the coefficients for all nonbasic variables in row i + // appear in Arow.row_start[i] to nonbasic_end[i]. + // Swap the coefficients for the leaving and entering variables + // The leaving variable is now nonbasic, the entering variable is now basic + row_mark.clear(); + const i_t col_start_enter = A.col_start[entering]; + const i_t col_end_enter = A.col_start[entering + 1]; + for (i_t p = col_start_enter; p < col_end_enter; p++) { + const i_t i = A.i[p]; + row_mark.push_back(i); + } + work_estimate += 2 * (col_end_enter - col_start_enter); + + // Move the coefficients for the entering variable to the end of the nonbasics + // And decrement the nonbasic count for the row + for (i_t i : row_mark) { + const i_t row_start = Arow.row_start[i]; + const i_t nb_end = nonbasic_end[i]; + for (i_t p = row_start; p <= nb_end; p++) { + const i_t j = Arow.j[p]; + if (j == entering) { + std::swap(Arow.j[p], Arow.j[nb_end]); + std::swap(Arow.x[p], Arow.x[nb_end]); + nonbasic_end[i]--; + break; + } + } + work_estimate += nb_end - row_start; + } + work_estimate += 2 * row_mark.size(); + + row_mark.clear(); + const i_t col_start_leave = A.col_start[leaving]; + const i_t col_end_leave = A.col_start[leaving + 1]; + for (i_t p = col_start_leave; p < col_end_leave; p++) { + const i_t i = A.i[p]; + row_mark.push_back(i); + } + work_estimate += 2 * (col_end_leave - col_start_leave); + + // Move the coefficient for the leaving variable to the end of the nonbasics + // And increment the nonbasic count for the row + for (i_t i : row_mark) { + const i_t nb_end = nonbasic_end[i]; + const i_t row_end = Arow.row_start[i + 1]; + for (i_t p = nb_end + 1; p < row_end; p++) { + const i_t j = Arow.j[p]; + if (j == leaving) { + std::swap(Arow.j[p], Arow.j[nb_end + 1]); + std::swap(Arow.x[p], Arow.x[nb_end + 1]); + nonbasic_end[i]++; + break; + } + } + work_estimate += row_end - nb_end; + } + work_estimate += 2 * row_mark.size(); } namespace phase2 { @@ -807,7 +864,7 @@ void update_single_primal_infeasibility(const std::vector& lower, } template -void update_primal_infeasibilities(const lp_problem_t& lp, +bool update_primal_infeasibilities(const lp_problem_t& lp, const simplex_solver_settings_t& settings, const std::vector& basic_list, const std::vector& x, @@ -819,6 +876,7 @@ void update_primal_infeasibilities(const lp_problem_t& lp, f_t& primal_inf, f_t& work_estimate) { + bool became_feasible = false; const f_t primal_tol = settings.primal_tol; const i_t nz = basic_change_list.size(); for (i_t k = 0; k < nz; ++k) { @@ -831,8 +889,10 @@ void update_primal_infeasibilities(const lp_problem_t& lp, const f_t old_val = squared_infeasibilities[j]; squared_infeasibilities[j] = 0.0; primal_inf = std::max(0.0, primal_inf - old_val); + if (old_val != 0.0) { became_feasible = true; } continue; } + const f_t old_val = squared_infeasibilities[j]; update_single_primal_infeasibility(lp.lower, lp.upper, x, @@ -841,8 +901,10 @@ void update_primal_infeasibilities(const lp_problem_t& lp, infeasibility_indices, j, primal_inf); + if (old_val != 0.0 && squared_infeasibilities[j] == 0.0) { became_feasible = true; } } work_estimate += 8 * nz; + return became_feasible; } template @@ -850,33 +912,14 @@ void clean_up_infeasibilities(std::vector& squared_infeasibilities, std::vector& infeasibility_indices, f_t& work_estimate) { - bool needs_clean_up = false; - const i_t initial_nz = infeasibility_indices.size(); - for (i_t k = 0; k < initial_nz; ++k) { - const i_t j = infeasibility_indices[k]; - const f_t squared_infeas = squared_infeasibilities[j]; - if (squared_infeas == 0.0) { needs_clean_up = true; } - } - work_estimate += 2 * initial_nz; - - if (needs_clean_up) { - i_t num_cleans = 0; - work_estimate += 2 * infeasibility_indices.size(); - for (size_t k = 0; k < infeasibility_indices.size(); ++k) { - const i_t j = infeasibility_indices[k]; - const f_t squared_infeas = squared_infeasibilities[j]; - if (squared_infeas == 0.0) { - const i_t new_j = infeasibility_indices.back(); - infeasibility_indices[k] = new_j; - infeasibility_indices.pop_back(); - if (squared_infeasibilities[new_j] == 0.0) { - k--; - } // Decrement k so that we process the same index again - num_cleans++; - } - } - work_estimate += 4 * num_cleans; + i_t write = 0; + const i_t n = infeasibility_indices.size(); + for (i_t k = 0; k < n; ++k) { + const i_t j = infeasibility_indices[k]; + if (squared_infeasibilities[j] != 0.0) { infeasibility_indices[write++] = j; } } + infeasibility_indices.resize(write); + work_estimate += 3 * n; } template @@ -1480,8 +1523,10 @@ i_t compute_steepest_edge_norm_entering(const simplex_solver_settings_t& ft, i_t basic_leaving_index, i_t entering_index, + i_t leaving_index, std::vector& steepest_edge_norms) { +#if 0 sparse_vector_t es_sparse(m, 1); es_sparse.i[0] = basic_leaving_index; es_sparse.x[0] = -1.0; @@ -1489,6 +1534,16 @@ i_t compute_steepest_edge_norm_entering(const simplex_solver_settings_t 1e-1) { + settings.log.printf("Steepest edge norm %e for entering %d. Leaving %d norm %e. Diff %e\n", + steepest_edge_norms[entering_index], + entering_index, + leaving_index, + steepest_edge_norms[leaving_index], + std::abs(steepest_edge_norms[entering_index] - steepest_edge_norms[leaving_index])); + } +#endif + steepest_edge_norms[entering_index] = steepest_edge_norms[leaving_index]; #ifdef STEEPEST_EDGE_DEBUG settings.log.printf("Steepest edge norm %e for entering j %d at i %d\n", steepest_edge_norms[entering_index], @@ -2408,6 +2463,7 @@ class phase2_timers_t { se_norms_time(0), se_entering_time(0), lu_update_time(0), + lu_factorization_time(0), perturb_time(0), vector_time(0), objective_time(0), @@ -2431,8 +2487,9 @@ class phase2_timers_t { { if (!record_time) { return; } const f_t total_time = bfrt_time + pricing_time + btran_time + ftran_time + flip_time + - delta_z_time + lu_update_time + se_norms_time + se_entering_time + - perturb_time + vector_time + objective_time + update_infeasibility_time; + delta_z_time + lu_update_time + lu_factorization_time + se_norms_time + + se_entering_time + perturb_time + vector_time + objective_time + + update_infeasibility_time; // clang-format off settings.log.printf("BFRT time %.2fs %4.1f%\n", bfrt_time, 100.0 * bfrt_time / total_time); settings.log.printf("Pricing time %.2fs %4.1f%\n", pricing_time, 100.0 * pricing_time / total_time); @@ -2441,6 +2498,7 @@ class phase2_timers_t { settings.log.printf("Flip time %.2fs %4.1f%\n", flip_time, 100.0 * flip_time / total_time); settings.log.printf("Delta_z time %.2fs %4.1f%\n", delta_z_time, 100.0 * delta_z_time / total_time); settings.log.printf("LU update time %.2fs %4.1f%\n", lu_update_time, 100.0 * lu_update_time / total_time); + settings.log.printf("LU factor time %.2fs %4.1f%\n", lu_factorization_time, 100.0 * lu_factorization_time / total_time); settings.log.printf("SE norms time %.2fs %4.1f%\n", se_norms_time, 100.0 * se_norms_time / total_time); settings.log.printf("SE enter time %.2fs %4.1f%\n", se_entering_time, 100.0 * se_entering_time / total_time); settings.log.printf("Perturb time %.2fs %4.1f%\n", perturb_time, 100.0 * perturb_time / total_time); @@ -2459,6 +2517,7 @@ class phase2_timers_t { f_t se_norms_time; f_t se_entering_time; f_t lu_update_time; + f_t lu_factorization_time; f_t perturb_time; f_t vector_time; f_t objective_time; @@ -2561,6 +2620,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, phase2::bound_info(lp, settings); phase2_work_estimate += 2 * n; + f_t refactor_work = 0.0; + f_t solve_work = 0.0; + if (initialize_basis) { PHASE2_NVTX_RANGE("DualSimplex::init_basis"); std::vector superbasic_list; @@ -2573,8 +2635,10 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, assert(superbasic_list.size() == 0); assert(nonbasic_list.size() == n - m); - i_t refactor_status = ft.refactor_basis( + f_t refactor_start_work = ft.work_estimate(); + i_t refactor_status = ft.refactor_basis( lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); + refactor_work = ft.work_estimate() - refactor_start_work; if (refactor_status == CONCURRENT_HALT_RETURN) { return dual::status_t::CONCURRENT_LIMIT; } if (refactor_status == TIME_LIMIT_RETURN) { return dual::status_t::TIME_LIMIT; } if (refactor_status > 0) { return dual::status_t::NUMERICAL; } @@ -2735,9 +2799,14 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 0); #endif - csc_matrix_t A_transpose(1, 1, 0); - lp.A.transpose(A_transpose); + csr_matrix_t Arow(1, 1, 0); + lp.A.to_compressed_row(Arow); phase2_work_estimate += 2 * lp.A.col_start[lp.A.n]; + std::vector nonbasic_end(m); + std::vector row_mark; + row_mark.reserve(m); + compute_initial_nonbasic_end(basic_mark, Arow, nonbasic_end); + phase2_work_estimate += lp.A.col_start[lp.A.n]; f_t obj = compute_objective(lp, x); phase2_work_estimate += 2 * n; @@ -2749,7 +2818,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, i_t num_refactors = 0; i_t total_bound_flips = 0; f_t delta_y_nz_percentage = 0.0; - phase2::phase2_timers_t timers(false); + phase2::phase2_timers_t timers(true); // Sparse vectors for main loop (declared outside loop for instrumentation) sparse_vector_t delta_y_sparse(m, 0); @@ -2781,6 +2850,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, 0.0, toc(start_time)); } + i_t iterations_since_refactor = 0; while (iter < iter_limit) { PHASE2_NVTX_RANGE("DualSimplex::phase2_main_loop"); @@ -2903,11 +2973,13 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, timers.start_timer(); delta_y_sparse.clear(); UTsol_sparse.clear(); + f_t btran_start_work = ft.work_estimate(); { PHASE2_NVTX_RANGE("DualSimplex::btran"); phase2::compute_delta_y(ft, basic_leaving_index, direction, delta_y_sparse, UTsol_sparse); } timers.btran_time += timers.stop_timer(); + solve_work += (ft.work_estimate() - btran_start_work); const f_t steepest_edge_norm_check = delta_y_sparse.norm2_squared(); phase2_work_estimate += 2 * delta_y_sparse.i.size(); @@ -2939,11 +3011,11 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, PHASE2_NVTX_RANGE("DualSimplex::delta_z"); if (use_transpose) { sparse_delta_z++; - compute_delta_z(A_transpose, + compute_delta_z(Arow, delta_y_sparse, leaving_index, direction, - nonbasic_mark, + nonbasic_end, delta_z_mark, delta_z_indices, delta_z, @@ -3278,6 +3350,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, utilde_sparse.clear(); scaled_delta_xB_sparse.clear(); rhs_sparse.from_csc_column(lp.A, entering_index); + f_t ftran_start_work = ft.work_estimate(); { PHASE2_NVTX_RANGE("DualSimplex::ftran"); if (phase2::compute_delta_x(lp, @@ -3299,7 +3372,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, return dual::status_t::NUMERICAL; } } - + solve_work += (ft.work_estimate() - ftran_start_work); timers.ftran_time += timers.stop_timer(); #ifdef CHECK_PRIMAL_STEP @@ -3310,6 +3383,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, #endif timers.start_timer(); + f_t se_norms_start_work = ft.work_estimate(); const i_t steepest_edge_status = phase2::update_steepest_edge_norms(settings, basic_list, ft, @@ -3331,6 +3405,7 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, #endif assert(steepest_edge_status == 0); timers.se_norms_time += timers.stop_timer(); + solve_work += (ft.work_estimate() - se_norms_start_work); timers.start_timer(); // x <- x + delta_x @@ -3365,42 +3440,48 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, #ifdef CHECK_BASIC_INFEASIBILITIES phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 2); #endif - phase2::update_primal_infeasibilities(lp, - settings, - basic_list, - x, - entering_index, - leaving_index, - delta_xB_0_sparse.i, - squared_infeasibilities, - infeasibility_indices, - primal_infeasibility_squared, - phase2_work_estimate); + bool needs_cleanup = phase2::update_primal_infeasibilities(lp, + settings, + basic_list, + x, + entering_index, + leaving_index, + delta_xB_0_sparse.i, + squared_infeasibilities, + infeasibility_indices, + primal_infeasibility_squared, + phase2_work_estimate); // Update primal infeasibilities due to changes in basic variables // from the leaving and entering variables - phase2::update_primal_infeasibilities(lp, - settings, - basic_list, - x, - entering_index, - leaving_index, - scaled_delta_xB_sparse.i, - squared_infeasibilities, - infeasibility_indices, - primal_infeasibility_squared, - phase2_work_estimate); + needs_cleanup |= phase2::update_primal_infeasibilities(lp, + settings, + basic_list, + x, + entering_index, + leaving_index, + scaled_delta_xB_sparse.i, + squared_infeasibilities, + infeasibility_indices, + primal_infeasibility_squared, + phase2_work_estimate); // Update the entering variable - phase2::update_single_primal_infeasibility(lp.lower, - lp.upper, - x, - settings.primal_tol, - squared_infeasibilities, - infeasibility_indices, - entering_index, - primal_infeasibility_squared); - - phase2::clean_up_infeasibilities( - squared_infeasibilities, infeasibility_indices, phase2_work_estimate); + { + const f_t old_val = squared_infeasibilities[entering_index]; + phase2::update_single_primal_infeasibility(lp.lower, + lp.upper, + x, + settings.primal_tol, + squared_infeasibilities, + infeasibility_indices, + entering_index, + primal_infeasibility_squared); + needs_cleanup |= (old_val != 0.0 && squared_infeasibilities[entering_index] == 0.0); + } + + if (needs_cleanup) { + phase2::clean_up_infeasibilities( + squared_infeasibilities, infeasibility_indices, phase2_work_estimate); + } #if CHECK_PRIMAL_INFEASIBILITIES phase2::check_primal_infeasibilities( @@ -3432,6 +3513,9 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, basic_mark[leaving_index] = -1; basic_mark[entering_index] = basic_leaving_index; + update_Arow( + leaving_index, entering_index, lp.A, row_mark, nonbasic_end, Arow, phase2_work_estimate); + #ifdef CHECK_BASIC_INFEASIBILITIES phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 5); #endif @@ -3440,22 +3524,35 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, // Refactor or update the basis factorization { PHASE2_NVTX_RANGE("DualSimplex::basis_update"); - bool should_refactor = ft.num_updates() > settings.refactor_frequency; + iterations_since_refactor++; + bool should_refactor = + (ft.num_updates() > 100 && solve_work > refactor_work) || (ft.num_updates() > 1000); + // settings.log.printf("Solve work %e refactor work %e iterations since refactor %d\n", + // solve_work, refactor_work, iterations_since_refactor); if (!should_refactor) { i_t recommend_refactor = ft.update(utilde_sparse, UTsol_sparse, basic_leaving_index); #ifdef CHECK_UPDATE phase2::check_update(lp, settings, ft, basic_list, basic_leaving_index); #endif should_refactor = recommend_refactor == 1; + timers.lu_update_time += timers.stop_timer(); + timers.start_timer(); } #ifdef CHECK_BASIC_INFEASIBILITIES phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 6); #endif if (should_refactor) { + settings.log.printf( + "Refactoring basis. Iteration %d solve work %e refactor work %e updates %d\n", + iter, + solve_work, + refactor_work, + ft.num_updates()); PHASE2_NVTX_RANGE("DualSimplex::refactorization"); num_refactors++; - bool should_recompute_x = true; // Need for numerically difficult problems like cbs-cta + bool should_recompute_x = true; // Needed for numerically difficult problems like cbs-cta + f_t refactor_start_work = ft.work_estimate(); i_t refactor_status = ft.refactor_basis( lp.A, settings, lp.lower, lp.upper, start_time, basic_list, nonbasic_list, vstatus); if (refactor_status == CONCURRENT_HALT_RETURN) { return dual::status_t::CONCURRENT_LIMIT; } @@ -3488,9 +3585,11 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, settings.log.printf("Successfully repaired basis. Iteration %d\n", iter); } + refactor_work = ft.work_estimate() - refactor_start_work; phase2::reset_basis_mark( basic_list, nonbasic_list, basic_mark, nonbasic_mark, phase2_work_estimate); + compute_initial_nonbasic_end(basic_mark, Arow, nonbasic_end); if (should_recompute_x) { std::vector unperturbed_x(n); phase2_work_estimate += n; @@ -3514,16 +3613,18 @@ dual::status_t dual_phase2_with_advanced_basis(i_t phase, infeasibility_indices, primal_infeasibility); phase2_work_estimate += 4 * m + 2 * n; + iterations_since_refactor = 0; + solve_work = 0.0; } #ifdef CHECK_BASIC_INFEASIBILITIES phase2::check_basic_infeasibilities(basic_list, basic_mark, infeasibility_indices, 7); #endif } - timers.lu_update_time += timers.stop_timer(); + timers.lu_factorization_time += timers.stop_timer(); timers.start_timer(); phase2::compute_steepest_edge_norm_entering( - settings, m, ft, basic_leaving_index, entering_index, delta_y_steepest_edge); + settings, m, ft, basic_leaving_index, entering_index, leaving_index, delta_y_steepest_edge); timers.se_entering_time += timers.stop_timer(); #ifdef STEEPEST_EDGE_DEBUG @@ -3653,15 +3754,19 @@ template void compute_reduced_cost_update(const lp_problem_t& delta_z, double& work_estimate); -template void compute_delta_z(const csc_matrix_t& A_transpose, +template void compute_delta_z(const csr_matrix_t& Arow, const sparse_vector_t& delta_y, int leaving_index, int direction, - const std::vector& nonbasic_mark, + const std::vector& nonbasic_end, std::vector& delta_z_mark, std::vector& delta_z_indices, std::vector& delta_z, double& work_estimate); + +template void compute_initial_nonbasic_end(const std::vector& basic_mark, + csr_matrix_t& Arow, + std::vector& nonbasic_end); #endif } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/dual_simplex/phase2.hpp b/cpp/src/dual_simplex/phase2.hpp index 5db797449..266d57c19 100644 --- a/cpp/src/dual_simplex/phase2.hpp +++ b/cpp/src/dual_simplex/phase2.hpp @@ -94,14 +94,19 @@ void compute_reduced_cost_update(const lp_problem_t& lp, f_t& work_estimate); template -void compute_delta_z(const csc_matrix_t& A_transpose, +void compute_delta_z(const csr_matrix_t& Arow, const sparse_vector_t& delta_y, i_t leaving_index, i_t direction, - const std::vector& nonbasic_mark, + const std::vector& nonbasic_end, std::vector& delta_z_mark, std::vector& delta_z_indices, std::vector& delta_z, f_t& work_estimate); +template +void compute_initial_nonbasic_end(const std::vector& basic_mark, + csr_matrix_t& Arow, + std::vector& nonbasic_end); + } // namespace cuopt::linear_programming::dual_simplex