From ff8213a6a9ff22839354c7f15740b911419460c7 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Tue, 2 Feb 2016 13:00:45 -0800 Subject: [PATCH] dev(lp): diminish the number of pivots of primal by sorting the columns Signed-off-by: Lev Nachmanson --- src/util/lp/lp_dual_core_solver.cpp | 8 ------ src/util/lp/lp_dual_core_solver.h | 2 -- src/util/lp/lp_primal_core_solver.cpp | 40 ++++++++++++++++++++++----- src/util/lp/lp_primal_core_solver.h | 4 ++- src/util/lp/lu.cpp | 2 +- 5 files changed, 37 insertions(+), 19 deletions(-) diff --git a/src/util/lp/lp_dual_core_solver.cpp b/src/util/lp/lp_dual_core_solver.cpp index 37c7ffb09..bf70d4d76 100644 --- a/src/util/lp/lp_dual_core_solver.cpp +++ b/src/util/lp/lp_dual_core_solver.cpp @@ -56,14 +56,6 @@ template void lp_dual_core_solver::fill_non_basis } } -template void lp_dual_core_solver::print_nb(std::ostream & out) { - out << "this is nb " << std::endl; - for (auto l : this->m_factorization->m_non_basic_columns) { - out << l << " "; - } - out << std::endl; -} - template void lp_dual_core_solver::restore_non_basis() { auto & nb = this->m_factorization->m_non_basic_columns; nb.clear(); diff --git a/src/util/lp/lp_dual_core_solver.h b/src/util/lp/lp_dual_core_solver.h index 87bbcfac0..0484d7f9d 100644 --- a/src/util/lp/lp_dual_core_solver.h +++ b/src/util/lp/lp_dual_core_solver.h @@ -49,8 +49,6 @@ public: void fill_non_basis_with_only_able_to_enter_columns(); - void print_nb(std::ostream & out); - void restore_non_basis(); bool update_basis(int entering, int leaving); diff --git a/src/util/lp/lp_primal_core_solver.cpp b/src/util/lp/lp_primal_core_solver.cpp index 3dfda91f3..0bf48e73a 100644 --- a/src/util/lp/lp_primal_core_solver.cpp +++ b/src/util/lp/lp_primal_core_solver.cpp @@ -8,17 +8,41 @@ #include "util/lp/lp_primal_core_solver.h" namespace lean { + // This core solver solves (Ax=b, low_bound_values \leq x \leq upper_bound_values, maximize costs*x ) // The right side b is given implicitly by x and the basis +template +void lp_primal_core_solver::sort_non_basis() { + for (unsigned j : this->m_non_basic_columns) { + T const & da = this->m_d[j]; + m_steepest_edge_coefficients[j] = da * da / this->m_column_norms[j]; + } + + std::sort(this->m_non_basic_columns.begin(), this->m_non_basic_columns.end(), [this](unsigned a, unsigned b) { + return m_steepest_edge_coefficients[a] > m_steepest_edge_coefficients[b]; + }); + // reinit m_basis_heading + for (int j = 0; j < this->m_non_basic_columns.size(); j++) + this->m_basis_heading[this->m_non_basic_columns[j]] = - j - 1; +} + + template int lp_primal_core_solver::choose_entering_column(unsigned number_of_benefitial_columns_to_go_over) { // at this moment m_y = cB * B(-1) lean_assert(number_of_benefitial_columns_to_go_over); - unsigned offset_in_nb = my_random() % this->m_non_basic_columns.size(); + if (m_sort_columns_counter == 0) { + sort_non_basis(); + m_sort_columns_counter = 20; + } else { + m_sort_columns_counter--; + } + lean_assert(m_sort_columns_counter>=0); + unsigned offset_in_nb = 0; unsigned initial_offset_in_non_basis = offset_in_nb; T steepest_edge = zero_of_type(); int entering = -1; do { - unsigned j = this->m_factorization->m_non_basic_columns[offset_in_nb]; + unsigned j = this->m_non_basic_columns[offset_in_nb]; push_forward_offset_in_non_basis(offset_in_nb); if (m_forbidden_enterings.find(j) != m_forbidden_enterings.end()) { continue; @@ -225,7 +249,8 @@ template lp_primal_core_solver:: lp_primal_core column_type_array, low_bound_values, upper_bound_values), - m_beta(A.row_count()) { + m_beta(A.row_count()), + m_steepest_edge_coefficients(A.column_count()) { this->m_status = UNKNOWN; this->m_column_norm_update_counter = settings.column_norms_update_frequency; } @@ -250,7 +275,8 @@ lp_primal_core_solver(static_matrix & A, column_type_array, m_low_bound_values_dummy, upper_bound_values), - m_beta(A.row_count()) { + m_beta(A.row_count()), + m_steepest_edge_coefficients(A.column_count()) { m_converted_harris_eps = convert_struct::convert(this->m_settings.harris_feasibility_tolerance); lean_assert(initial_x_is_correct()); m_low_bound_values_dummy.resize(A.column_count(), zero_of_type()); @@ -504,8 +530,8 @@ template void lp_primal_core_solver::push_forw offset_in_nb = 0; } -template unsigned lp_primal_core_solver::get_number_of_non_basic_column_to_try_for_enter() { - unsigned ret = this->m_factorization->m_non_basic_columns.size(); +template unsigned lp_primal_core_solver::get_number_of_non_basic_column_to_try_for_enter() { + unsigned ret = this->m_non_basic_columns.size(); if (this->m_status == TENTATIVE_UNBOUNDED) return ret; // we really need to find entering with a large reduced cost if (ret > 300) { @@ -594,7 +620,7 @@ template void lp_primal_core_solver::delete_fa } // according to Swietanowski, " A new steepest edge approximation for the simplex method for linear programming" -template void lp_primal_core_solver::init_column_norms() { +template void lp_primal_core_solver::init_column_norms() { m_column_norm_update_counter = 0; for (unsigned j : this->m_non_basic_columns) this->m_column_norms[j] = 1; diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 4629e54d1..1528e731a 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -48,7 +48,9 @@ public: std::vector m_costs_backup; bool m_current_x_is_feasible; T m_converted_harris_eps = convert_struct::convert(this->m_settings.harris_feasibility_tolerance); - + unsigned m_sort_columns_counter = 0; + std::vector m_steepest_edge_coefficients; + void sort_non_basis(); int choose_entering_column(unsigned number_of_benefitial_columns_to_go_over); int find_leaving_and_t_precisely(unsigned entering, X & t); diff --git a/src/util/lp/lu.cpp b/src/util/lp/lu.cpp index 85ba4375b..6761b9d5f 100644 --- a/src/util/lp/lu.cpp +++ b/src/util/lp/lu.cpp @@ -646,7 +646,7 @@ row_eta_matrix *lu::get_row_eta_matrix_and_set_row_vector(unsigned r #endif !m_settings.abs_val_is_smaller_than_pivot_tolerance((m_row_eta_work_vector[lowest_row_of_the_bump] - pivot_elem_for_checking) / denom)) { set_status(LU_status::Degenerated); - std::cout << "diagonal element is off" << std::endl; + // std::cout << "diagonal element is off" << std::endl; return nullptr; } #ifdef LEAN_DEBUG