diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt index 5719de44fae..cc58c9610fc 100644 --- a/src/math/lp/CMakeLists.txt +++ b/src/math/lp/CMakeLists.txt @@ -19,7 +19,6 @@ z3_add_component(lp lar_solver.cpp lar_core_solver.cpp lp_core_solver_base.cpp - lp_dual_core_solver.cpp lp_primal_core_solver.cpp lp_settings.cpp lu.cpp diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index 508a06d730e..0d9f4297419 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -27,7 +27,6 @@ template bool lp::lp_core_solver_base::basis_heading_is_correct( template void lp::lp_core_solver_base::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned); template bool lp::lp_core_solver_base::column_is_dual_feasible(unsigned int) const; template void lp::lp_core_solver_base::fill_reduced_costs_from_m_y_by_rows(); -template bool lp::lp_core_solver_base::find_x_by_solving(); template lp::non_basic_column_value_position lp::lp_core_solver_base::get_non_basic_column_value_position(unsigned int) const; template lp::non_basic_column_value_position lp::lp_core_solver_base >::get_non_basic_column_value_position(unsigned int) const; template lp::non_basic_column_value_position lp::lp_core_solver_base::get_non_basic_column_value_position(unsigned int) const; @@ -53,7 +52,6 @@ template bool lp::lp_core_solver_base::basis_heading_is_correc template void lp::lp_core_solver_base::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned); template bool lp::lp_core_solver_base::column_is_dual_feasible(unsigned int) const; template void lp::lp_core_solver_base::fill_reduced_costs_from_m_y_by_rows(); -template bool lp::lp_core_solver_base::find_x_by_solving(); template bool lp::lp_core_solver_base::print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(char const*, std::ostream &); template void lp::lp_core_solver_base::restore_x(unsigned int, lp::mpq const&); template void lp::lp_core_solver_base::set_non_basic_x_to_correct_bounds(); @@ -112,7 +110,6 @@ template bool lp::lp_core_solver_base::column_is_feasible(unsi template bool lp::lp_core_solver_base >::column_is_feasible(unsigned int) const; template bool lp::lp_core_solver_base >::snap_non_basic_x_to_bound(); template void lp::lp_core_solver_base >::init_lu(); -template bool lp::lp_core_solver_base >::find_x_by_solving(); template void lp::lp_core_solver_base >::restore_x(unsigned int, lp::numeric_pair const&); template bool lp::lp_core_solver_base>::pivot_column_tableau(unsigned int, unsigned int); template bool lp::lp_core_solver_base::pivot_column_tableau(unsigned int, unsigned int); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index edff349d23d..729792c1266 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -298,9 +298,6 @@ class lp_core_solver_base { void rs_minus_Anx(vector & rs); - bool find_x_by_solving(); - - bool basis_has_no_doubles() const; bool non_basis_has_no_doubles() const; diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index ae04bd5eec6..e470072e04c 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -338,12 +338,6 @@ rs_minus_Anx(vector & rs) { } } -template bool lp_core_solver_base:: -find_x_by_solving() { - solve_Ax_eq_b(); - return true; -} - template bool lp_core_solver_base::column_is_feasible(unsigned j) const { const X& x = this->m_x[j]; switch (this->m_column_types[j]) { diff --git a/src/math/lp/lp_dual_core_solver.cpp b/src/math/lp/lp_dual_core_solver.cpp deleted file mode 100644 index 8cf45f7c6cd..00000000000 --- a/src/math/lp/lp_dual_core_solver.cpp +++ /dev/null @@ -1,44 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#include -#include -#include -#include "util/vector.h" -#include -#include "math/lp/lp_dual_core_solver_def.h" -template void lp::lp_dual_core_solver::start_with_initial_basis_and_make_it_dual_feasible(); -template void lp::lp_dual_core_solver::solve(); -template lp::lp_dual_core_solver::lp_dual_core_solver(lp::static_matrix&, vector&, - vector&, - vector&, - vector&, - vector &, - vector &, - vector&, - vector&, - vector&, - vector&, - lp::lp_settings&, const lp::column_namer&); -template void lp::lp_dual_core_solver::start_with_initial_basis_and_make_it_dual_feasible(); -template void lp::lp_dual_core_solver::solve(); -template void lp::lp_dual_core_solver::restore_non_basis(); -template void lp::lp_dual_core_solver::restore_non_basis(); -template void lp::lp_dual_core_solver::revert_to_previous_basis(); -template void lp::lp_dual_core_solver::revert_to_previous_basis(); diff --git a/src/math/lp/lp_dual_core_solver.h b/src/math/lp/lp_dual_core_solver.h deleted file mode 100644 index b78a2214701..00000000000 --- a/src/math/lp/lp_dual_core_solver.h +++ /dev/null @@ -1,212 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once -#include "math/lp/static_matrix.h" -#include "math/lp/lp_core_solver_base.h" -#include -#include -#include -#include -#include "util/vector.h" - -namespace lp { -template -class lp_dual_core_solver:public lp_core_solver_base { -public: - vector & m_can_enter_basis; - int m_r; // the row of the leaving column - int m_p; // leaving column; that is m_p = m_basis[m_r] - T m_delta; // the offset of the leaving basis variable - int m_sign_of_alpha_r; // see page 27 - T m_theta_D; - T m_theta_P; - int m_q; - // todo : replace by a vector later - std::set m_breakpoint_set; // it is F in "Progress in the dual simplex method ..." - std::set m_flipped_boxed; - std::set m_tight_set; // it is the set of all breakpoints that become tight when m_q becomes tight - vector m_a_wave; - vector m_betas; // m_betas[i] is approximately a square of the norm of the i-th row of the reverse of B - T m_harris_tolerance; - std::set m_forbidden_rows; - - lp_dual_core_solver(static_matrix & A, - vector & can_enter_basis, - vector & b, // the right side vector - vector & x, // the number of elements in x needs to be at least as large as the number of columns in A - vector & basis, - vector & nbasis, - vector & heading, - vector & costs, - vector & column_type_array, - vector & lower_bound_values, - vector & upper_bound_values, - lp_settings & settings, - const column_namer & column_names): - lp_core_solver_base(A, - // b, - basis, - nbasis, - heading, - x, - costs, - settings, - column_names, - column_type_array, - lower_bound_values, - upper_bound_values), - m_can_enter_basis(can_enter_basis), - m_a_wave(this->m_m()), - m_betas(this->m_m()) { - m_harris_tolerance = numeric_traits::precise()? numeric_traits::zero() : T(this->m_settings.harris_feasibility_tolerance); - lp_assert(false); - this->init_basic_part_of_basis_heading(); - fill_non_basis_with_only_able_to_enter_columns(); - } - - void init_a_wave_by_zeros(); - - void fill_non_basis_with_only_able_to_enter_columns() { - auto & nb = this->m_nbasis; - nb.reset(); - unsigned j = this->m_n(); - while (j--) { - if (this->m_basis_heading[j] >= 0 || !m_can_enter_basis[j]) continue; - nb.push_back(j); - this->m_basis_heading[j] = - static_cast(nb.size()); - } - } - - void restore_non_basis(); - - bool update_basis(int entering, int leaving); - - void recalculate_xB_and_d(); - - void recalculate_d(); - - void init_betas(); - - void adjust_xb_for_changed_xn_and_init_betas(); - - void start_with_initial_basis_and_make_it_dual_feasible(); - - bool done(); - - T get_edge_steepness_for_lower_bound(unsigned p); - - T get_edge_steepness_for_upper_bound(unsigned p); - - T pricing_for_row(unsigned i); - - void pricing_loop(unsigned number_of_rows_to_try, unsigned offset_in_rows); - - bool advance_on_known_p(); - - int define_sign_of_alpha_r(); - - bool can_be_breakpoint(unsigned j); - - void fill_breakpoint_set(); - - void DSE_FTran(); - T get_delta(); - - void restore_d(); - - bool d_is_correct(); - - void xb_minus_delta_p_pivot_column(); - - void update_betas(); - - void apply_flips(); - - void snap_xN_column_to_bounds(unsigned j); - - void snap_xN_to_bounds(); - - void init_beta_precisely(unsigned i); - - void init_betas_precisely(); - - // step 7 of the algorithm from Progress - bool basis_change_and_update(); - - void revert_to_previous_basis(); - - non_basic_column_value_position m_entering_boundary_position; - bool update_basis_and_x_local(int entering, int leaving, X const & tt); - void recover_leaving(); - - bool problem_is_dual_feasible() const; - - bool snap_runaway_nonbasic_column(unsigned); - - bool snap_runaway_nonbasic_columns(); - - unsigned get_number_of_rows_to_try_for_leaving(); - - void update_a_wave(const T & del, unsigned j) { - this->m_A.add_column_to_vector(del, j, & m_a_wave[0]); - } - - bool delta_keeps_the_sign(int initial_delta_sign, const T & delta); - - void set_status_to_tentative_dual_unbounded_or_dual_unbounded(); - - // it is positive if going from low bound to upper bound and negative if going from upper bound to low bound - T signed_span_of_boxed(unsigned j) { - return this->x_is_at_lower_bound(j)? this->bound_span(j): - this->bound_span(j); - } - - void add_tight_breakpoints_and_q_to_flipped_set(); - - T delta_lost_on_flips_of_tight_breakpoints(); - - bool tight_breakpoinst_are_all_boxed(); - - T calculate_harris_delta_on_breakpoint_set(); - - void fill_tight_set_on_harris_delta(const T & harris_delta ); - - void find_q_on_tight_set(); - - void find_q_and_tight_set(); - - void erase_tight_breakpoints_and_q_from_breakpoint_set(); - - bool ratio_test(); - - void process_flipped(); - void update_d_and_xB(); - - void calculate_beta_r_precisely(); - // see "Progress in the dual simplex method for large scale LP problems: practical dual phase 1 algorithms" - - void update_xb_after_bound_flips(); - - void one_iteration(); - - void solve(); - - bool lower_bounds_are_set() const override { return true; } -}; -} diff --git a/src/math/lp/lp_dual_core_solver_def.h b/src/math/lp/lp_dual_core_solver_def.h deleted file mode 100644 index 9c39fe195bf..00000000000 --- a/src/math/lp/lp_dual_core_solver_def.h +++ /dev/null @@ -1,699 +0,0 @@ -/*++ -Copyright (c) 2017 Microsoft Corporation - -Module Name: - - - -Abstract: - - - -Author: - - Lev Nachmanson (levnach) - -Revision History: - - ---*/ -#pragma once - -#include -#include -#include "util/vector.h" -#include "math/lp/lp_dual_core_solver.h" - -namespace lp { - -template void lp_dual_core_solver::init_a_wave_by_zeros() { - unsigned j = this->m_m(); - while (j--) { - m_a_wave[j] = numeric_traits::zero(); - } -} - -template void lp_dual_core_solver::restore_non_basis() { - auto & nb = this->m_nbasis; - nb.reset(); - unsigned j = this->m_n(); - while (j--) { - if (this->m_basis_heading[j] >= 0 ) continue; - if (m_can_enter_basis[j]) { - lp_assert(std::find(nb.begin(), nb.end(), j) == nb.end()); - nb.push_back(j); - this->m_basis_heading[j] = - static_cast(nb.size()); - } - } -} - -template bool lp_dual_core_solver::update_basis(int entering, int leaving) { - // the second argument is the element of the entering column from the pivot row - its value should be equal to the low diagonal element of the bump after all pivoting is done - if (this->m_refactor_counter++ < 200) { - this->m_factorization->replace_column(this->m_ed[this->m_factorization->basis_heading(leaving)], this->m_w); - if (this->m_factorization->get_status() == LU_status::OK) { - this->m_factorization->change_basis(entering, leaving); - return true; - } - } - // need to refactor - this->m_factorization->change_basis(entering, leaving); - init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings); - this->m_refactor_counter = 0; - if (this->m_factorization->get_status() != LU_status::OK) { - LP_OUT(this->m_settings, "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->total_iterations() << std::endl); - this->m_iters_with_no_cost_growing++; - return false; - } - return true; -} - -template void lp_dual_core_solver::recalculate_xB_and_d() { - this->solve_Ax_eq_b(); - recalculate_d(); -} - -template void lp_dual_core_solver::recalculate_d() { -lp_assert(false) -} - -template void lp_dual_core_solver::init_betas() { - // todo : look at page 194 of Progress in the dual simplex algorithm for solving large scale LP problems : techniques for a fast and stable implementation - // the current implementation is not good enough: todo - unsigned i = this->m_m(); - while (i--) { - m_betas[i] = 1; - } -} - -template void lp_dual_core_solver::adjust_xb_for_changed_xn_and_init_betas() { - this->solve_Ax_eq_b(); - init_betas(); -} - -template void lp_dual_core_solver::start_with_initial_basis_and_make_it_dual_feasible() { - this->set_non_basic_x_to_correct_bounds(); // It is not an efficient version, see 3.29, - // however this version does not require that m_x is the solution of Ax = 0 beforehand - adjust_xb_for_changed_xn_and_init_betas(); -} - -template bool lp_dual_core_solver::done() { - if (this->get_status() == lp_status::OPTIMAL) { - return true; - } - - return false; // todo, need to be more cases -} - -template T lp_dual_core_solver::get_edge_steepness_for_lower_bound(unsigned p) { - lp_assert(this->m_basis_heading[p] >= 0 && static_cast(this->m_basis_heading[p]) < this->m_m()); - T del = this->m_x[p] - this->m_lower_bounds[p]; - del *= del; - return del / this->m_betas[this->m_basis_heading[p]]; -} - -template T lp_dual_core_solver::get_edge_steepness_for_upper_bound(unsigned p) { - lp_assert(this->m_basis_heading[p] >= 0 && static_cast(this->m_basis_heading[p]) < this->m_m()); - T del = this->m_x[p] - this->m_upper_bounds[p]; - del *= del; - return del / this->m_betas[this->m_basis_heading[p]]; -} - -template T lp_dual_core_solver::pricing_for_row(unsigned i) { - unsigned p = this->m_basis[i]; - switch (this->m_column_types[p]) { - case column_type::fixed: - case column_type::boxed: - if (this->x_below_low_bound(p)) { - T del = get_edge_steepness_for_lower_bound(p); - return del; - } - if (this->x_above_upper_bound(p)) { - T del = get_edge_steepness_for_upper_bound(p); - return del; - } - return numeric_traits::zero(); - case column_type::lower_bound: - if (this->x_below_low_bound(p)) { - T del = get_edge_steepness_for_lower_bound(p); - return del; - } - return numeric_traits::zero(); - break; - case column_type::upper_bound: - if (this->x_above_upper_bound(p)) { - T del = get_edge_steepness_for_upper_bound(p); - return del; - } - return numeric_traits::zero(); - break; - case column_type::free_column: - lp_assert(numeric_traits::is_zero(this->m_d[p])); - return numeric_traits::zero(); - default: - lp_unreachable(); - } - lp_unreachable(); - return numeric_traits::zero(); -} - -template void lp_dual_core_solver::pricing_loop(unsigned number_of_rows_to_try, unsigned offset_in_rows) { - m_r = -1; - T steepest_edge_max = numeric_traits::zero(); - unsigned initial_offset_in_rows = offset_in_rows; - unsigned i = offset_in_rows; - unsigned rows_left = number_of_rows_to_try; - do { - if (m_forbidden_rows.find(i) != m_forbidden_rows.end()) { - if (++i == this->m_m()) { - i = 0; - } - continue; - } - T se = pricing_for_row(i); - if (se > steepest_edge_max) { - steepest_edge_max = se; - m_r = i; - if (rows_left > 0) { - rows_left--; - } - } - if (++i == this->m_m()) { - i = 0; - } - } while (i != initial_offset_in_rows && rows_left); - if (m_r == -1) { - if (this->get_status() != lp_status::UNSTABLE) { - this->set_status(lp_status::OPTIMAL); - } - } else { - m_p = this->m_basis[m_r]; - m_delta = get_delta(); - if (advance_on_known_p()){ - m_forbidden_rows.clear(); - return; - } - // failure in advance_on_known_p - if (this->get_status() == lp_status::FLOATING_POINT_ERROR) { - return; - } - this->set_status(lp_status::UNSTABLE); - m_forbidden_rows.insert(m_r); - } -} - - // this calculation is needed for the steepest edge update, - // it hijackes m_pivot_row_of_B_1 for this purpose since we will need it anymore to the end of the cycle -template void lp_dual_core_solver::DSE_FTran() { // todo, see algorithm 7 from page 35 - this->m_factorization->solve_By_for_T_indexed_only(this->m_pivot_row_of_B_1, this->m_settings); -} - -template bool lp_dual_core_solver::advance_on_known_p() { - - return false; -} - -template int lp_dual_core_solver::define_sign_of_alpha_r() { - switch (this->m_column_types[m_p]) { - case column_type::boxed: - case column_type::fixed: - if (this->x_below_low_bound(m_p)) { - return -1; - } - if (this->x_above_upper_bound(m_p)) { - return 1; - } - lp_unreachable(); - case column_type::lower_bound: - if (this->x_below_low_bound(m_p)) { - return -1; - } - lp_unreachable(); - case column_type::upper_bound: - if (this->x_above_upper_bound(m_p)) { - return 1; - } - lp_unreachable(); - default: - lp_unreachable(); - } - lp_unreachable(); - return 0; -} - -template bool lp_dual_core_solver::can_be_breakpoint(unsigned j) { - if (this->pivot_row_element_is_too_small_for_ratio_test(j)) return false; - switch (this->m_column_types[j]) { - case column_type::lower_bound: - lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_lower_bounds[j])); - return m_sign_of_alpha_r * this->m_pivot_row[j] > 0; - case column_type::upper_bound: - lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_upper_bounds[j])); - return m_sign_of_alpha_r * this->m_pivot_row[j] < 0; - case column_type::boxed: - { - bool lower_bound = this->x_is_at_lower_bound(j); - bool grawing = m_sign_of_alpha_r * this->m_pivot_row[j] > 0; - return lower_bound == grawing; - } - case column_type::fixed: // is always dual feasible so we ignore it - return false; - case column_type::free_column: - return true; - default: - return false; - } -} - -template void lp_dual_core_solver::fill_breakpoint_set() { - m_breakpoint_set.clear(); - for (unsigned j : this->non_basis()) { - if (can_be_breakpoint(j)) { - m_breakpoint_set.insert(j); - } - } -} - -// template void lp_dual_core_solver::FTran() { -// this->solve_Bd(m_q); -// } - -template T lp_dual_core_solver::get_delta() { - switch (this->m_column_types[m_p]) { - case column_type::boxed: - if (this->x_below_low_bound(m_p)) { - return this->m_x[m_p] - this->m_lower_bounds[m_p]; - } - if (this->x_above_upper_bound(m_p)) { - return this->m_x[m_p] - this->m_upper_bounds[m_p]; - } - lp_unreachable(); - case column_type::lower_bound: - if (this->x_below_low_bound(m_p)) { - return this->m_x[m_p] - this->m_lower_bounds[m_p]; - } - lp_unreachable(); - case column_type::upper_bound: - if (this->x_above_upper_bound(m_p)) { - return get_edge_steepness_for_upper_bound(m_p); - } - lp_unreachable(); - case column_type::fixed: - return this->m_x[m_p] - this->m_upper_bounds[m_p]; - default: - lp_unreachable(); - } - lp_unreachable(); - return zero_of_type(); -} - -template void lp_dual_core_solver::restore_d() { - this->m_d[m_p] = numeric_traits::zero(); - for (auto j : this->non_basis()) { - this->m_d[j] += m_theta_D * this->m_pivot_row[j]; - } -} - -template bool lp_dual_core_solver::d_is_correct() { - - lp_assert(false); - return true; -} - -template void lp_dual_core_solver::xb_minus_delta_p_pivot_column() { - unsigned i = this->m_m(); - while (i--) { - this->m_x[this->m_basis[i]] -= m_theta_P * this->m_ed[i]; - } -} - -template void lp_dual_core_solver::update_betas() { // page 194 of Progress ... todo - once in a while betas have to be reinitialized - T one_over_arq = numeric_traits::one() / this->m_pivot_row[m_q]; - T beta_r = this->m_betas[m_r] = std::max(T(0.0001), (m_betas[m_r] * one_over_arq) * one_over_arq); - T k = -2 * one_over_arq; - unsigned i = this->m_m(); - while (i--) { - if (static_cast(i) == m_r) continue; - T a = this->m_ed[i]; - m_betas[i] += a * (a * beta_r + k * this->m_pivot_row_of_B_1[i]); - if (m_betas[i] < T(0.0001)) - m_betas[i] = T(0.0001); - } -} - -template void lp_dual_core_solver::apply_flips() { - for (unsigned j : m_flipped_boxed) { - lp_assert(this->x_is_at_bound(j)); - if (this->x_is_at_lower_bound(j)) { - this->m_x[j] = this->m_upper_bounds[j]; - } else { - this->m_x[j] = this->m_lower_bounds[j]; - } - } -} - -template void lp_dual_core_solver::snap_xN_column_to_bounds(unsigned j) { - switch (this->m_column_type[j]) { - case column_type::fixed: - this->m_x[j] = this->m_lower_bounds[j]; - break; - case column_type::boxed: - if (this->x_is_at_lower_bound(j)) { - this->m_x[j] = this->m_lower_bounds[j]; - } else { - this->m_x[j] = this->m_upper_bounds[j]; - } - break; - case column_type::lower_bound: - this->m_x[j] = this->m_lower_bounds[j]; - break; - case column_type::upper_bound: - this->m_x[j] = this->m_upper_bounds[j]; - break; - case column_type::free_column: - break; - default: - lp_unreachable(); - } -} - -template void lp_dual_core_solver::snap_xN_to_bounds() { - for (auto j : this->non_basis()) { - snap_xN_column_to_bounds(j); - } -} - -template void lp_dual_core_solver::init_beta_precisely(unsigned i) { - vector vec(this->m_m(), numeric_traits::zero()); - vec[i] = numeric_traits::one(); - this->m_factorization->solve_yB_with_error_check(vec, this->m_basis); - T beta = numeric_traits::zero(); - for (T & v : vec) { - beta += v * v; - } - this->m_betas[i] =beta; -} - -template void lp_dual_core_solver::init_betas_precisely() { - unsigned i = this->m_m(); - while (i--) { - init_beta_precisely(i); - } -} - -// step 7 of the algorithm from Progress -template bool lp_dual_core_solver::basis_change_and_update() { - - return true; -} - -template void lp_dual_core_solver::recover_leaving() { - switch (m_entering_boundary_position) { - case at_lower_bound: - case at_fixed: - this->m_x[m_q] = this->m_lower_bounds[m_q]; - break; - case at_upper_bound: - this->m_x[m_q] = this->m_upper_bounds[m_q]; - break; - case free_of_bounds: - this->m_x[m_q] = zero_of_type(); - default: - lp_unreachable(); - } -} - -template void lp_dual_core_solver::revert_to_previous_basis() { - LP_OUT(this->m_settings, "revert to previous basis on ( " << m_p << ", " << m_q << ")" << std::endl); - this->change_basis_unconditionally(m_p, m_q); - init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_settings); - if (this->m_factorization->get_status() != LU_status::OK) { - this->set_status(lp_status::FLOATING_POINT_ERROR); // complete failure - return; - } - recover_leaving(); - if (!this->find_x_by_solving()) { - this->set_status(lp_status::FLOATING_POINT_ERROR); - return; - } - recalculate_xB_and_d(); - init_betas_precisely(); -} - -// returns true if the column has been snapped -template bool lp_dual_core_solver::snap_runaway_nonbasic_column(unsigned j) { - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::lower_bound: - if (!this->x_is_at_lower_bound(j)) { - this->m_x[j] = this->m_lower_bounds[j]; - return true; - } - break; - case column_type::boxed: - { - bool closer_to_lower_bound = abs(this->m_lower_bounds[j] - this->m_x[j]) < abs(this->m_upper_bounds[j] - this->m_x[j]); - if (closer_to_lower_bound) { - if (!this->x_is_at_lower_bound(j)) { - this->m_x[j] = this->m_lower_bounds[j]; - return true; - } - } else { - if (!this->x_is_at_upper_bound(j)) { - this->m_x[j] = this->m_lower_bounds[j]; - return true; - } - } - } - break; - case column_type::upper_bound: - if (!this->x_is_at_upper_bound(j)) { - this->m_x[j] = this->m_upper_bounds[j]; - return true; - } - break; - default: - break; - } - return false; -} - - -template bool lp_dual_core_solver::problem_is_dual_feasible() const { - for (unsigned j : this->non_basis()){ - if (!this->column_is_dual_feasible(j)) { - return false; - } - } - return true; -} - -template unsigned lp_dual_core_solver::get_number_of_rows_to_try_for_leaving() { - unsigned s = this->m_m(); - if (this->m_m() > 300) { - s = (unsigned)((s / 100.0) * this->m_settings.percent_of_entering_to_check); - } - return this->m_settings.random_next() % s + 1; -} - -template bool lp_dual_core_solver::delta_keeps_the_sign(int initial_delta_sign, const T & delta) { - if (numeric_traits::precise()) - return ((delta > numeric_traits::zero()) && (initial_delta_sign == 1)) || - ((delta < numeric_traits::zero()) && (initial_delta_sign == -1)); - - double del = numeric_traits::get_double(delta); - return ( (del > this->m_settings.zero_tolerance) && (initial_delta_sign == 1)) || - ((del < - this->m_settings.zero_tolerance) && (initial_delta_sign == -1)); -} - -template void lp_dual_core_solver::set_status_to_tentative_dual_unbounded_or_dual_unbounded() { - if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) { - this->set_status(lp_status::DUAL_UNBOUNDED); - } else { - this->set_status(lp_status::TENTATIVE_DUAL_UNBOUNDED); - } -} - -template void lp_dual_core_solver::add_tight_breakpoints_and_q_to_flipped_set() { - m_flipped_boxed.insert(m_q); - for (auto j : m_tight_set) { - m_flipped_boxed.insert(j); - } -} - -template T lp_dual_core_solver::delta_lost_on_flips_of_tight_breakpoints() { - T ret = abs(this->bound_span(m_q) * this->m_pivot_row[m_q]); - for (auto j : m_tight_set) { - ret += abs(this->bound_span(j) * this->m_pivot_row[j]); - } - return ret; -} - -template bool lp_dual_core_solver::tight_breakpoinst_are_all_boxed() { - if (this->m_column_types[m_q] != column_type::boxed) return false; - for (auto j : m_tight_set) { - if (this->m_column_types[j] != column_type::boxed) return false; - } - return true; -} - -template T lp_dual_core_solver::calculate_harris_delta_on_breakpoint_set() { - bool first_time = true; - T ret = zero_of_type(); - lp_assert(m_breakpoint_set.size() > 0); - for (auto j : m_breakpoint_set) { - T t; - if (this->x_is_at_lower_bound(j)) { - t = abs((std::max(this->m_d[j], numeric_traits::zero()) + m_harris_tolerance) / this->m_pivot_row[j]); - } else { - t = abs((std::min(this->m_d[j], numeric_traits::zero()) - m_harris_tolerance) / this->m_pivot_row[j]); - } - if (first_time) { - ret = t; - first_time = false; - } else if (t < ret) { - ret = t; - } - } - return ret; -} - -template void lp_dual_core_solver::fill_tight_set_on_harris_delta(const T & harris_delta ){ - m_tight_set.clear(); - for (auto j : m_breakpoint_set) { - if (this->x_is_at_lower_bound(j)) { - if (abs(std::max(this->m_d[j], numeric_traits::zero()) / this->m_pivot_row[j]) <= harris_delta){ - m_tight_set.insert(j); - } - } else { - if (abs(std::min(this->m_d[j], numeric_traits::zero() ) / this->m_pivot_row[j]) <= harris_delta){ - m_tight_set.insert(j); - } - } - } -} - -template void lp_dual_core_solver::find_q_on_tight_set() { - m_q = -1; - T max_pivot; - for (auto j : m_tight_set) { - T r = abs(this->m_pivot_row[j]); - if (m_q != -1) { - if (r > max_pivot) { - max_pivot = r; - m_q = j; - } - } else { - max_pivot = r; - m_q = j; - } - } - m_tight_set.erase(m_q); - lp_assert(m_q != -1); -} - -template void lp_dual_core_solver::find_q_and_tight_set() { - T harris_del = calculate_harris_delta_on_breakpoint_set(); - fill_tight_set_on_harris_delta(harris_del); - find_q_on_tight_set(); - m_entering_boundary_position = this->get_non_basic_column_value_position(m_q); -} - -template void lp_dual_core_solver::erase_tight_breakpoints_and_q_from_breakpoint_set() { - m_breakpoint_set.erase(m_q); - for (auto j : m_tight_set) { - m_breakpoint_set.erase(j); - } -} - -template bool lp_dual_core_solver::ratio_test() { - m_sign_of_alpha_r = define_sign_of_alpha_r(); - fill_breakpoint_set(); - m_flipped_boxed.clear(); - int initial_delta_sign = m_delta >= numeric_traits::zero()? 1: -1; - do { - if (m_breakpoint_set.empty()) { - set_status_to_tentative_dual_unbounded_or_dual_unbounded(); - return false; - } - this->set_status(lp_status::FEASIBLE); - find_q_and_tight_set(); - if (!tight_breakpoinst_are_all_boxed()) break; - T del = m_delta - delta_lost_on_flips_of_tight_breakpoints() * initial_delta_sign; - if (!delta_keeps_the_sign(initial_delta_sign, del)) break; - if (m_tight_set.size() + 1 == m_breakpoint_set.size()) { - break; // deciding not to flip since we might get stuck without finding m_q, the column entering the basis - } - // we can flip m_q together with the tight set and look for another breakpoint candidate for m_q and another tight set - add_tight_breakpoints_and_q_to_flipped_set(); - m_delta = del; - erase_tight_breakpoints_and_q_from_breakpoint_set(); - } while (true); - m_theta_D = this->m_d[m_q] / this->m_pivot_row[m_q]; - return true; -} - -template void lp_dual_core_solver::process_flipped() { - init_a_wave_by_zeros(); - for (auto j : m_flipped_boxed) { - update_a_wave(signed_span_of_boxed(j), j); - } -} -template void lp_dual_core_solver::update_d_and_xB() { - for (auto j : this->non_basis()) { - this->m_d[j] -= m_theta_D * this->m_pivot_row[j]; - } - this->m_d[m_p] = - m_theta_D; - if (!m_flipped_boxed.empty()) { - process_flipped(); - update_xb_after_bound_flips(); - } -} - -template void lp_dual_core_solver::calculate_beta_r_precisely() { - T t = numeric_traits::zero(); - unsigned i = this->m_m(); - while (i--) { - T b = this->m_pivot_row_of_B_1[i]; - t += b * b; - } - m_betas[m_r] = t; -} -// see "Progress in the dual simplex method for large scale LP problems: practical dual phase 1 algorithms" - -template void lp_dual_core_solver::update_xb_after_bound_flips() { - this->m_factorization->solve_By(m_a_wave); - unsigned i = this->m_m(); - while (i--) { - this->m_x[this->m_basis[i]] -= m_a_wave[i]; - } -} - -template void lp_dual_core_solver::one_iteration() { - unsigned number_of_rows_to_try = get_number_of_rows_to_try_for_leaving(); - unsigned offset_in_rows = this->m_settings.random_next() % this->m_m(); - if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) { - number_of_rows_to_try = this->m_m(); - } else { - this->set_status(lp_status::FEASIBLE); - } - pricing_loop(number_of_rows_to_try, offset_in_rows); - lp_assert(problem_is_dual_feasible()); -} - -template void lp_dual_core_solver::solve() { // see the page 35 - lp_assert(d_is_correct()); - lp_assert(problem_is_dual_feasible()); - lp_assert(this->basis_heading_is_correct()); - //this->set_total_iterations(0); - this->iters_with_no_cost_growing() = 0; - do { - if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over("", *this->m_settings.get_message_ostream())){ - return; - } - one_iteration(); - } while (this->get_status() != lp_status::FLOATING_POINT_ERROR && this->get_status() != lp_status::DUAL_UNBOUNDED && this->get_status() != lp_status::OPTIMAL && - this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements - ); -} -} diff --git a/src/math/lp/static_matrix.cpp b/src/math/lp/static_matrix.cpp index cde96277b5e..93360b9140f 100644 --- a/src/math/lp/static_matrix.cpp +++ b/src/math/lp/static_matrix.cpp @@ -23,7 +23,6 @@ Revision History: #include #include "math/lp/static_matrix_def.h" #include "math/lp/lp_core_solver_base.h" -#include "math/lp/lp_dual_core_solver.h" #include "math/lp/lp_primal_core_solver.h" #include "math/lp/scaler.h" #include "math/lp/lar_solver.h"