diff --git a/src/ast/rewriter/seq_rewriter.cpp b/src/ast/rewriter/seq_rewriter.cpp index 6fd518cb0..24df51845 100644 --- a/src/ast/rewriter/seq_rewriter.cpp +++ b/src/ast/rewriter/seq_rewriter.cpp @@ -5024,12 +5024,14 @@ br_status seq_rewriter::mk_re_star(expr* a, expr_ref& result) { * (re.range c_1 c_n) */ br_status seq_rewriter::mk_re_range(expr* lo, expr* hi, expr_ref& result) { - zstring s; + zstring slo, shi; unsigned len = 0; bool is_empty = false; - if (str().is_string(lo, s) && s.length() != 1) + if (str().is_string(lo, slo) && slo.length() != 1) is_empty = true; - if (str().is_string(hi, s) && s.length() != 1) + if (str().is_string(hi, shi) && shi.length() != 1) + is_empty = true; + if (slo.length() == 1 && shi.length() == 1 && slo[0] > shi[0]) is_empty = true; len = min_length(lo).second; if (len > 1) @@ -5246,7 +5248,17 @@ br_status seq_rewriter::reduce_re_is_empty(expr* r, expr_ref& result) { else if (re().is_range(r, r1, r2) && str().is_string(r1, s1) && str().is_string(r2, s2) && s1.length() == 1 && s2.length() == 1) { - result = m().mk_bool_val(s1[0] <= s2[0]); + result = m().mk_bool_val(s1[0] > s2[0]); + return BR_DONE; + } + else if (re().is_range(r, r1, r2) && + str().is_string(r1, s1) && s1.length() != 1) { + result = m().mk_true(); + return BR_DONE; + } + else if (re().is_range(r, r1, r2) && + str().is_string(r2, s2) && s2.length() != 1) { + result = m().mk_true(); return BR_DONE; } else if ((re().is_loop(r, r1, lo) || @@ -5307,6 +5319,7 @@ br_status seq_rewriter::mk_le_core(expr * l, expr * r, expr_ref & result) { } br_status seq_rewriter::mk_eq_core(expr * l, expr * r, expr_ref & result) { + TRACE("seq", tout << mk_pp(l, m()) << " == " << mk_pp(r, m()) << "\n"); expr_ref_vector res(m()); expr_ref_pair_vector new_eqs(m()); if (m_util.is_re(l)) { diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp index 69a9a101c..48015471c 100644 --- a/src/math/lp/int_solver.cpp +++ b/src/math/lp/int_solver.cpp @@ -14,34 +14,51 @@ namespace lp { -int_solver::patcher::patcher(int_solver& lia): - lia(lia), - lra(lia.lra), - lrac(lia.lrac) -{} - - -lia_move int_solver::patcher::patch_nbasic_columns() { - lia.settings().stats().m_patches++; - lp_assert(lia.is_feasible()); - m_patch_success = 0; - m_patch_fail = 0; - unsigned num = lra.A_r().column_count(); - for (unsigned v = 0; v < num; v++) { - if (lia.is_base(v)) - continue; - patch_nbasic_column(v); + int_solver::patcher::patcher(int_solver& lia): + lia(lia), + lra(lia.lra), + lrac(lia.lrac) + {} + + void int_solver::patcher::remove_fixed_vars_from_base() { + unsigned num = lra.A_r().column_count(); + for (unsigned v = 0; v < num; v++) { + if (!lia.is_base(v) || !lia.is_fixed(v)) + continue; + auto const & r = lra.basic2row(v); + for (auto const& c : r) { + if (c.var() != v && !lia.is_fixed(c.var())) { + lra.pivot(c.var(), v); + break; + } + } + } } -// for (unsigned j : lia.lrac.m_r_nbasis) -// patch_nbasic_column(j); - lp_assert(lia.is_feasible()); - verbose_stream() << "patch " << m_patch_success << " fails " << m_patch_fail << "\n"; - if (!lia.has_inf_int()) { - lia.settings().stats().m_patches_success++; - return lia_move::sat; + + lia_move int_solver::patcher::patch_nbasic_columns() { + remove_fixed_vars_from_base(); + lia.settings().stats().m_patches++; + lp_assert(lia.is_feasible()); + m_patch_success = 0; + m_patch_fail = 0; + m_num_ones = 0; + m_num_divides = 0; + unsigned num = lra.A_r().column_count(); + for (unsigned v = 0; v < num; v++) { + if (lia.is_base(v)) + continue; + patch_nbasic_column(v); + } + lp_assert(lia.is_feasible()); + verbose_stream() << "patch " << m_patch_success << " fails " << m_patch_fail << " ones " << m_num_ones << " divides " << m_num_divides << "\n"; + //lra.display(verbose_stream()); + //exit(0); + if (!lia.has_inf_int()) { + lia.settings().stats().m_patches_success++; + return lia_move::sat; + } + return lia_move::undef; } - return lia_move::undef; -} void int_solver::patcher::patch_nbasic_column(unsigned j) { impq & val = lrac.m_r_x[j]; @@ -54,8 +71,24 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) { bool m_is_one = m.is_one(); bool val_is_int = lia.value_is_int(j); + + const auto & A = lra.A_r(); + if (!m_is_one) { + verbose_stream() << "divisor: " << m << "\n"; + lia.display_column(verbose_stream(), j); + for (auto c : A.column(j)) { + unsigned row_index = c.var(); + lia.display_row(verbose_stream(), A.m_rows[row_index]); + verbose_stream() << "\n"; + } + } + // check whether value of j is already a multiple of m. if (val_is_int && (m_is_one || (val.x / m).is_int())) { + if (m_is_one) + ++m_num_ones; + else + ++m_num_divides; return; } @@ -77,13 +110,15 @@ void int_solver::patcher::patch_nbasic_column(unsigned j) { } #endif +#if 0 verbose_stream() << "TARGET v" << j << " -> ["; if (inf_l) verbose_stream() << "-oo"; else verbose_stream() << ceil(l.x); verbose_stream() << ", "; if (inf_u) verbose_stream() << "oo"; else verbose_stream() << floor(u.x); verbose_stream() << "]"; verbose_stream() << ", m: " << m << ", val: " << val << ", is_int: " << lra.column_is_int(j) << "\n"; - +#endif + if (!inf_l) { l = impq(m_is_one ? ceil(l) : m * ceil(l / m)); if (inf_u || l <= u) { @@ -472,14 +507,11 @@ bool int_solver::at_upper(unsigned j) const { } std::ostream & int_solver::display_row(std::ostream & out, lp::row_strip const & row) const { -bool first = true; + bool first = true; auto & rslv = lrac.m_r_solver; -for (const auto &c : row) - { - if (is_fixed(c.var())) - { - if (!get_value(c.var()).is_zero()) - { + for (const auto &c : row) { + if (is_fixed(c.var())) { + if (!get_value(c.var()).is_zero()) { impq val = get_value(c.var()) * c.coeff(); if (!first && val.is_pos()) out << "+"; @@ -498,17 +530,11 @@ for (const auto &c : row) } else if (c.coeff().is_minus_one()) out << "-"; - else - { - if (c.coeff().is_pos()) - { - if (!first) - out << "+"; - } + else { + if (c.coeff().is_pos() && !first) + out << "+"; if (c.coeff().is_big()) - { out << " b*"; - } else out << c.coeff(); } @@ -516,8 +542,7 @@ for (const auto &c : row) first = false; } out << "\n"; - for (const auto &c : row) - { + for (const auto &c : row) { if (is_fixed(c.var())) continue; rslv.print_column_info(c.var(), out); @@ -526,14 +551,13 @@ for (const auto &c : row) } return out; } + std::ostream& int_solver::display_row_info(std::ostream & out, unsigned row_index) const { auto & rslv = lrac.m_r_solver; auto const& row = rslv.m_A.m_rows[row_index]; return display_row(out, row); } - - bool int_solver::shift_var(unsigned j, unsigned range) { if (is_fixed(j) || is_base(j)) return false; diff --git a/src/math/lp/int_solver.h b/src/math/lp/int_solver.h index 7f1807058..9a7d75573 100644 --- a/src/math/lp/int_solver.h +++ b/src/math/lp/int_solver.h @@ -46,12 +46,15 @@ class int_solver { lar_core_solver& lrac; unsigned m_patch_success = 0; unsigned m_patch_fail = 0; + unsigned m_num_ones = 0; + unsigned m_num_divides = 0; public: patcher(int_solver& lia); bool should_apply() const { return true; } lia_move operator()() { return patch_nbasic_columns(); } void patch_nbasic_column(unsigned j); private: + void remove_fixed_vars_from_base(); lia_move patch_nbasic_columns(); }; diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index 06ef4d50b..10cd53141 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -93,6 +93,8 @@ public: void solve(); + void pivot(int entering, int leaving) { m_r_solver.pivot(entering, leaving); } + bool lower_bounds_are_set() const { return true; } const indexed_vector & get_pivot_row() const { diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h index 762bae481..77fcbfc14 100644 --- a/src/math/lp/lar_solver.h +++ b/src/math/lp/lar_solver.h @@ -393,27 +393,32 @@ public: void clear_inf_set() { m_mpq_lar_core_solver.m_r_solver.inf_set().clear(); } + inline void remove_column_from_inf_set(unsigned j) { m_mpq_lar_core_solver.m_r_solver.remove_column_from_inf_set(j); } + + void pivot(int entering, int leaving) { + m_mpq_lar_core_solver.pivot(entering, leaving); + } + template void change_basic_columns_dependend_on_a_given_nb_column_report(unsigned j, const numeric_pair & delta, const ChangeReport& after) { - for (const auto & c : A_r().m_columns[j]) { - unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()]; - if (tableau_with_costs()) { - m_basic_columns_with_changed_cost.insert(bj); - } - m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, - A_r().get_val(c) * delta); - after(bj); - TRACE("change_x_del", - tout << "changed basis column " << bj << ", it is " << - ( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)? "feas":"inf") << std::endl;); - } - } - + for (const auto & c : A_r().m_columns[j]) { + unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()]; + if (tableau_with_costs()) + m_basic_columns_with_changed_cost.insert(bj); + m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, - A_r().get_val(c) * delta); + after(bj); + TRACE("change_x_del", + tout << "changed basis column " << bj << ", it is " << + ( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)?"feas":"inf") << std::endl;); + } + } + template void set_value_for_nbasic_column_report(unsigned j, const impq & new_val, diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index 587ac39a9..762beb289 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -324,7 +324,7 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { if(! m_A.pivot_row_to_row_given_cell(piv_row_index, c, j)) { return false; } - if (m_pivoted_rows!= nullptr) + if (m_pivoted_rows != nullptr) m_pivoted_rows->insert(c.var()); } diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h index 207428985..02a45bfd2 100644 --- a/src/math/lp/lp_primal_core_solver.h +++ b/src/math/lp/lp_primal_core_solver.h @@ -35,640 +35,624 @@ Revision History: #include namespace lp { -// This core solver solves (Ax=b, lower_bound_values \leq x \leq -// upper_bound_values, maximize costs*x ) The right side b is given implicitly -// by x and the basis -template -class lp_primal_core_solver : public lp_core_solver_base { -public: - int m_sign_of_entering_delta; - vector m_costs_backup; - unsigned m_inf_row_index_for_tableau; - bool m_bland_mode_tableau; - u_set m_left_basis_tableau; - unsigned m_bland_mode_threshold; - unsigned m_left_basis_repeated; - vector m_leaving_candidates; + // This core solver solves (Ax=b, lower_bound_values \leq x \leq + // upper_bound_values, maximize costs*x ) The right side b is given implicitly + // by x and the basis + template + class lp_primal_core_solver : public lp_core_solver_base { + public: + int m_sign_of_entering_delta; + vector m_costs_backup; + unsigned m_inf_row_index_for_tableau; + bool m_bland_mode_tableau; + u_set m_left_basis_tableau; + unsigned m_bland_mode_threshold; + unsigned m_left_basis_repeated; + vector m_leaving_candidates; - std::list m_non_basis_list; - void sort_non_basis(); - int choose_entering_column_tableau(); + std::list m_non_basis_list; + void sort_non_basis(); + int choose_entering_column_tableau(); - bool needs_to_grow(unsigned bj) const { - lp_assert(!this->column_is_feasible(bj)); - switch (this->m_column_types[bj]) { - case column_type::free_column: - return false; - case column_type::fixed: - case column_type::lower_bound: - case column_type::boxed: - return this->x_below_low_bound(bj); - default: - return false; - } - UNREACHABLE(); // unreachable - return false; - } + bool needs_to_grow(unsigned bj) const { + lp_assert(!this->column_is_feasible(bj)); + switch (this->m_column_types[bj]) { + case column_type::free_column: + return false; + case column_type::fixed: + case column_type::lower_bound: + case column_type::boxed: + return this->x_below_low_bound(bj); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } - int inf_sign_of_column(unsigned bj) const { - lp_assert(!this->column_is_feasible(bj)); - switch (this->m_column_types[bj]) { - case column_type::free_column: - return 0; - case column_type::lower_bound: - return 1; - case column_type::fixed: - case column_type::boxed: - return this->x_above_upper_bound(bj) ? -1 : 1; - default: - return -1; - } - UNREACHABLE(); // unreachable - return 0; - } + int inf_sign_of_column(unsigned bj) const { + lp_assert(!this->column_is_feasible(bj)); + switch (this->m_column_types[bj]) { + case column_type::free_column: + return 0; + case column_type::lower_bound: + return 1; + case column_type::fixed: + case column_type::boxed: + return this->x_above_upper_bound(bj) ? -1 : 1; + default: + return -1; + } + UNREACHABLE(); // unreachable + return 0; + } - bool monoid_can_decrease(const row_cell &rc) const { - unsigned j = rc.var(); - lp_assert(this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::free_column: - return true; - case column_type::fixed: - return false; - case column_type::lower_bound: - if (is_pos(rc.coeff())) { - return this->x_above_lower_bound(j); - } + bool monoid_can_decrease(const row_cell &rc) const { + unsigned j = rc.var(); + lp_assert(this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::free_column: + return true; + case column_type::fixed: + return false; + case column_type::lower_bound: + return !is_pos(rc.coeff()) || this->x_above_lower_bound(j); + case column_type::upper_bound: + return is_pos(rc.coeff()) || this->x_below_upper_bound(j); + case column_type::boxed: + if (is_pos(rc.coeff())) + return this->x_above_lower_bound(j); + return this->x_below_upper_bound(j); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } - return true; - case column_type::upper_bound: - if (is_pos(rc.coeff())) { - return true; - } + bool monoid_can_increase(const row_cell &rc) const { + unsigned j = rc.var(); + lp_assert(this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::free_column: + return true; + case column_type::fixed: + return false; + case column_type::lower_bound: + if (is_neg(rc.coeff())) + return this->x_above_lower_bound(j); + return true; + case column_type::upper_bound: + if (is_neg(rc.coeff())) + return true; + return this->x_below_upper_bound(j); + case column_type::boxed: + if (is_neg(rc.coeff())) + return this->x_above_lower_bound(j); + return this->x_below_upper_bound(j); + default: + return false; + } + UNREACHABLE(); // unreachable + return false; + } - return this->x_below_upper_bound(j); - case column_type::boxed: - if (is_pos(rc.coeff())) { - return this->x_above_lower_bound(j); - } + unsigned get_number_of_basic_vars_that_might_become_inf(unsigned j) const { + // consider looking at the signs here: todo + unsigned r = 0; + for (const auto &cc : this->m_A.m_columns[j]) { + unsigned k = this->m_basis[cc.var()]; + if (this->m_column_types[k] != column_type::free_column) + r++; + } + return r; + } - return this->x_below_upper_bound(j); - default: - return false; - } - UNREACHABLE(); // unreachable - return false; - } + int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T &a_ent) { + int j = -1; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (const row_cell &rc : this->m_A.m_rows[i]) { + if (rc.var() == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } + else { + if (!monoid_can_increase(rc)) + continue; + } + if (rc.var() < static_cast(j)) { + j = rc.var(); + a_ent = rc.coeff(); + } + } + if (j == -1) + m_inf_row_index_for_tableau = i; + return j; + } - bool monoid_can_increase(const row_cell &rc) const { - unsigned j = rc.var(); - lp_assert(this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::free_column: - return true; - case column_type::fixed: - return false; - case column_type::lower_bound: - if (is_neg(rc.coeff())) { - return this->x_above_lower_bound(j); - } + int find_beneficial_column_in_row_tableau_rows(int i, T &a_ent) { + if (m_bland_mode_tableau) + return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); + // a short row produces short infeasibility explanation and benefits at + // least one pivot operation + int choice = -1; + int nchoices = 0; + unsigned num_of_non_free_basics = 1000000; + unsigned len = 100000000; + unsigned bj = this->m_basis[i]; + bool bj_needs_to_grow = needs_to_grow(bj); + for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { + const row_cell &rc = this->m_A.m_rows[i][k]; + unsigned j = rc.var(); + if (j == bj) + continue; + if (bj_needs_to_grow) { + if (!monoid_can_decrease(rc)) + continue; + } else { + if (!monoid_can_increase(rc)) + continue; + } + unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); + if (damage < num_of_non_free_basics) { + num_of_non_free_basics = damage; + len = this->m_A.m_columns[j].size(); + choice = k; + nchoices = 1; + } + else if (damage == num_of_non_free_basics && + this->m_A.m_columns[j].size() <= len && + (this->m_settings.random_next() % (++nchoices))) { + choice = k; + len = this->m_A.m_columns[j].size(); + } + } - return true; - case column_type::upper_bound: - if (is_neg(rc.coeff())) { - return true; - } + if (choice == -1) { + m_inf_row_index_for_tableau = i; + return -1; + } + const row_cell &rc = this->m_A.m_rows[i][choice]; + a_ent = rc.coeff(); + return rc.var(); + } - return this->x_below_upper_bound(j); - case column_type::boxed: - if (is_neg(rc.coeff())) { - return this->x_above_lower_bound(j); - } + bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta, X &t, bool &unlimited); + + bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X &t); + + int find_leaving_and_t_tableau(unsigned entering, X &t); - return this->x_below_upper_bound(j); - default: - return false; - } - UNREACHABLE(); // unreachable - return false; - } + void limit_theta(const X &lim, X &theta, bool &unlimited) { + if (unlimited) { + theta = lim; + unlimited = false; + } + else + theta = std::min(lim, theta); + } - unsigned get_number_of_basic_vars_that_might_become_inf( - unsigned j) const { // consider looking at the signs here: todo - unsigned r = 0; - for (const auto &cc : this->m_A.m_columns[j]) { - unsigned k = this->m_basis[cc.var()]; - if (this->m_column_types[k] != column_type::free_column) - r++; - } - return r; - } + void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound); + limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited); + } - int find_beneficial_column_in_row_tableau_rows_bland_mode(int i, T &a_ent) { - int j = -1; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (const row_cell &rc : this->m_A.m_rows[i]) { - if (rc.var() == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) - continue; - } - if (rc.var() < static_cast(j)) { - j = rc.var(); - a_ent = rc.coeff(); - } - } - if (j == -1) { - m_inf_row_index_for_tableau = i; - } + void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound); + limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited); + } - return j; - } + void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound); + limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j], theta, unlimited); + } - int find_beneficial_column_in_row_tableau_rows(int i, T &a_ent) { - if (m_bland_mode_tableau) - return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); - // a short row produces short infeasibility explanation and benefits at - // least one pivot operation - int choice = -1; - int nchoices = 0; - unsigned num_of_non_free_basics = 1000000; - unsigned len = 100000000; - unsigned bj = this->m_basis[i]; - bool bj_needs_to_grow = needs_to_grow(bj); - for (unsigned k = 0; k < this->m_A.m_rows[i].size(); k++) { - const row_cell &rc = this->m_A.m_rows[i][k]; - unsigned j = rc.var(); - if (j == bj) - continue; - if (bj_needs_to_grow) { - if (!monoid_can_decrease(rc)) - continue; - } else { - if (!monoid_can_increase(rc)) - continue; - } - unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); - if (damage < num_of_non_free_basics) { - num_of_non_free_basics = damage; - len = this->m_A.m_columns[j].size(); - choice = k; - nchoices = 1; - } else if (damage == num_of_non_free_basics && - this->m_A.m_columns[j].size() <= len && - (this->m_settings.random_next() % (++nchoices))) { - choice = k; - len = this->m_A.m_columns[j].size(); - } - } + void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound); + limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, unlimited); + }; - if (choice == -1) { - m_inf_row_index_for_tableau = i; - return -1; - } - const row_cell &rc = this->m_A.m_rows[i][choice]; - a_ent = rc.coeff(); - return rc.var(); - } + void get_bound_on_variable_and_update_leaving_precisely( + unsigned j, vector &leavings, T m, X &t, + T &abs_of_d_of_leaving); - bool try_jump_to_another_bound_on_entering(unsigned entering, const X &theta, - X &t, bool &unlimited); - bool try_jump_to_another_bound_on_entering_unlimited(unsigned entering, X &t); - int find_leaving_and_t_tableau(unsigned entering, X &t); - - void limit_theta(const X &lim, X &theta, bool &unlimited) { - if (unlimited) { - theta = lim; - unlimited = false; - } else { - theta = std::min(lim, theta); - } - } - - void limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m < 0 && this->m_column_types[j] == column_type::upper_bound); - limit_inf_on_upper_bound_m_neg(m, this->m_x[j], this->m_upper_bounds[j], - theta, unlimited); - } - - void limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m < 0 && this->m_column_types[j] == column_type::lower_bound); - limit_inf_on_bound_m_neg(m, this->m_x[j], this->m_lower_bounds[j], theta, - unlimited); - } - - void limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m > 0 && this->m_column_types[j] == column_type::lower_bound); - limit_inf_on_lower_bound_m_pos(m, this->m_x[j], this->m_lower_bounds[j], - theta, unlimited); - } - - void limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m > 0 && this->m_column_types[j] == column_type::upper_bound); - limit_inf_on_bound_m_pos(m, this->m_x[j], this->m_upper_bounds[j], theta, - unlimited); - }; - - void get_bound_on_variable_and_update_leaving_precisely( - unsigned j, vector &leavings, T m, X &t, - T &abs_of_d_of_leaving); - - X get_max_bound(vector &b); + X get_max_bound(vector &b); #ifdef Z3DEBUG - void check_Ax_equal_b(); - void check_the_bounds(); - void check_bound(unsigned i); - void check_correctness(); + void check_Ax_equal_b(); + void check_the_bounds(); + void check_bound(unsigned i); + void check_correctness(); #endif - // from page 183 of Istvan Maros's book - // the basis structures have not changed yet - void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving); + // from page 183 of Istvan Maros's book + // the basis structures have not changed yet + void update_reduced_costs_from_pivot_row(unsigned entering, unsigned leaving); - // return 0 if the reduced cost at entering is close enough to the refreshed - // 1 if it is way off, and 2 if it is unprofitable - int refresh_reduced_cost_at_entering_and_check_that_it_is_off( - unsigned entering); + // return 0 if the reduced cost at entering is close enough to the refreshed + // 1 if it is way off, and 2 if it is unprofitable + int refresh_reduced_cost_at_entering_and_check_that_it_is_off( + unsigned entering); - void backup_and_normalize_costs(); + void backup_and_normalize_costs(); - void advance_on_entering_and_leaving_tableau(int entering, int leaving, X &t); - void advance_on_entering_equal_leaving_tableau(int entering, X &t); + void advance_on_entering_and_leaving_tableau(int entering, int leaving, X &t); + void advance_on_entering_equal_leaving_tableau(int entering, X &t); - bool need_to_switch_costs() const { - if (this->m_settings.simplex_strategy() == - simplex_strategy_enum::tableau_rows) - return false; - // lp_assert(calc_current_x_is_feasible() == - // current_x_is_feasible()); - return this->current_x_is_feasible() == this->using_infeas_costs(); - } - - void advance_on_entering_tableau(int entering); - - void push_forward_offset_in_non_basis(unsigned &offset_in_nb); - - unsigned get_number_of_non_basic_column_to_try_for_enter(); - - // returns the number of iterations - unsigned solve(); - - void find_feasible_solution(); - - // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;} - - void one_iteration_tableau(); - - // this version assumes that the leaving already has the right value, and does - // not update it - void update_x_tableau_rows(unsigned entering, unsigned leaving, - const X &delta) { - this->add_delta_to_x(entering, delta); - for (const auto &c : this->m_A.m_columns[entering]) { - if (leaving != this->m_basis[c.var()]) { - this->add_delta_to_x_and_track_feasibility( - this->m_basis[c.var()], -delta * this->m_A.get_val(c)); - } - } - } - - void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) { - lp_assert(entering != leaving); - update_x_tableau_rows(entering, leaving, tt); - this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); - this->change_basis(entering, leaving); - } - - void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, - const X &theta) { - update_basis_and_x_tableau_rows(entering, leaving, theta); - this->track_column_feasibility(entering); - } - - int find_smallest_inf_column() { - int j = -1; - for (unsigned k : this->inf_set()) { - if (k < static_cast(j)) { - j = k; - } - } - return j; - } - - const X &get_val_for_leaving(unsigned j) const { - lp_assert(!this->column_is_feasible(j)); - switch (this->m_column_types[j]) { - case column_type::fixed: - case column_type::upper_bound: - return this->m_upper_bounds[j]; - case column_type::lower_bound: - return this->m_lower_bounds[j]; - break; - case column_type::boxed: - if (this->x_above_upper_bound(j)) - return this->m_upper_bounds[j]; - else - return this->m_lower_bounds[j]; - break; - default: - UNREACHABLE(); - return this->m_lower_bounds[j]; - } - } - - void one_iteration_tableau_rows() { - int leaving = find_smallest_inf_column(); - if (leaving == -1) { - this->set_status(lp_status::OPTIMAL); - return; - } - - SASSERT(this->column_is_base(leaving)); - - if (!m_bland_mode_tableau) { - if (m_left_basis_tableau.contains(leaving)) { - if (++m_left_basis_repeated > m_bland_mode_threshold) { - m_bland_mode_tableau = true; + void pivot(int entering, int leaving) { + this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); + this->change_basis(entering, leaving); } - } else { - m_left_basis_tableau.insert(leaving); - } - } - T a_ent; - int entering = find_beneficial_column_in_row_tableau_rows( - this->m_basis_heading[leaving], a_ent); - if (entering == -1) { - this->set_status(lp_status::INFEASIBLE); - return; - } - const X &new_val_for_leaving = get_val_for_leaving(leaving); - X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; - this->m_x[leaving] = new_val_for_leaving; - this->remove_column_from_inf_set(leaving); - advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta); - if (this->current_x_is_feasible()) - this->set_status(lp_status::OPTIMAL); - } - void decide_on_status_when_cannot_find_entering() { - this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL - : lp_status::INFEASIBLE); - } + bool need_to_switch_costs() const { + if (this->m_settings.simplex_strategy() == + simplex_strategy_enum::tableau_rows) + return false; + // lp_assert(calc_current_x_is_feasible() == + // current_x_is_feasible()); + return this->current_x_is_feasible() == this->using_infeas_costs(); + } - void limit_theta_on_basis_column_for_feas_case_m_neg_no_check( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m < 0); - limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited); - if (theta < zero_of_type()) - theta = zero_of_type(); - } + void advance_on_entering_tableau(int entering); - bool limit_inf_on_bound_m_neg(const T &m, const X &x, const X &bound, - X &theta, bool &unlimited) { - // x gets smaller - lp_assert(m < 0); - if (this->below_bound(x, bound)) - return false; - if (this->above_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; - } - return true; - } + void push_forward_offset_in_non_basis(unsigned &offset_in_nb); - bool limit_inf_on_bound_m_pos(const T &m, const X &x, const X &bound, - X &theta, bool &unlimited) { - // x gets larger - lp_assert(m > 0); - if (this->above_bound(x, bound)) - return false; - if (this->below_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } else { - theta = zero_of_type(); - unlimited = false; - } + unsigned get_number_of_non_basic_column_to_try_for_enter(); - return true; - } + // returns the number of iterations + unsigned solve(); - void limit_inf_on_lower_bound_m_pos(const T &m, const X &x, const X &bound, + void find_feasible_solution(); + + // bool is_tiny() const {return this->m_m < 10 && this->m_n < 20;} + + void one_iteration_tableau(); + + // this version assumes that the leaving already has the right value, and does + // not update it + void update_x_tableau_rows(unsigned entering, unsigned leaving, + const X &delta) { + this->add_delta_to_x(entering, delta); + for (const auto &c : this->m_A.m_columns[entering]) + if (leaving != this->m_basis[c.var()]) + this->add_delta_to_x_and_track_feasibility( + this->m_basis[c.var()], -delta * this->m_A.get_val(c)); + } + + void update_basis_and_x_tableau_rows(int entering, int leaving, X const &tt) { + lp_assert(entering != leaving); + update_x_tableau_rows(entering, leaving, tt); + this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); + this->change_basis(entering, leaving); + } + + void advance_on_entering_and_leaving_tableau_rows(int entering, int leaving, + const X &theta) { + update_basis_and_x_tableau_rows(entering, leaving, theta); + this->track_column_feasibility(entering); + } + + int find_smallest_inf_column() { + int j = -1; + for (unsigned k : this->inf_set()) { + if (k < static_cast(j)) { + j = k; + } + } + return j; + } + + const X &get_val_for_leaving(unsigned j) const { + lp_assert(!this->column_is_feasible(j)); + switch (this->m_column_types[j]) { + case column_type::fixed: + case column_type::upper_bound: + return this->m_upper_bounds[j]; + case column_type::lower_bound: + return this->m_lower_bounds[j]; + break; + case column_type::boxed: + if (this->x_above_upper_bound(j)) + return this->m_upper_bounds[j]; + else + return this->m_lower_bounds[j]; + break; + default: + UNREACHABLE(); + return this->m_lower_bounds[j]; + } + } + + void one_iteration_tableau_rows() { + int leaving = find_smallest_inf_column(); + if (leaving == -1) { + this->set_status(lp_status::OPTIMAL); + return; + } + + SASSERT(this->column_is_base(leaving)); + + if (!m_bland_mode_tableau) { + if (m_left_basis_tableau.contains(leaving)) { + if (++m_left_basis_repeated > m_bland_mode_threshold) { + m_bland_mode_tableau = true; + } + } else { + m_left_basis_tableau.insert(leaving); + } + } + T a_ent; + int entering = find_beneficial_column_in_row_tableau_rows( + this->m_basis_heading[leaving], a_ent); + if (entering == -1) { + this->set_status(lp_status::INFEASIBLE); + return; + } + const X &new_val_for_leaving = get_val_for_leaving(leaving); + X theta = (this->m_x[leaving] - new_val_for_leaving) / a_ent; + this->m_x[leaving] = new_val_for_leaving; + this->remove_column_from_inf_set(leaving); + advance_on_entering_and_leaving_tableau_rows(entering, leaving, theta); + if (this->current_x_is_feasible()) + this->set_status(lp_status::OPTIMAL); + } + + void decide_on_status_when_cannot_find_entering() { + this->set_status(this->current_x_is_feasible() ? lp_status::OPTIMAL + : lp_status::INFEASIBLE); + } + + void limit_theta_on_basis_column_for_feas_case_m_neg_no_check( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m < 0); + limit_theta((this->m_lower_bounds[j] - this->m_x[j]) / m, theta, unlimited); + if (theta < zero_of_type()) + theta = zero_of_type(); + } + + bool limit_inf_on_bound_m_neg(const T &m, const X &x, const X &bound, X &theta, bool &unlimited) { - // x gets larger - lp_assert(m > 0); - if (this->below_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } - } + // x gets smaller + lp_assert(m < 0); + if (this->below_bound(x, bound)) + return false; + if (this->above_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } else { + theta = zero_of_type(); + unlimited = false; + } + return true; + } - void limit_inf_on_upper_bound_m_neg(const T &m, const X &x, const X &bound, + bool limit_inf_on_bound_m_pos(const T &m, const X &x, const X &bound, X &theta, bool &unlimited) { - // x gets smaller - lp_assert(m < 0); - if (this->above_bound(x, bound)) { - limit_theta((bound - x) / m, theta, unlimited); - } - } + // x gets larger + lp_assert(m > 0); + if (this->above_bound(x, bound)) + return false; + if (this->below_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } else { + theta = zero_of_type(); + unlimited = false; + } - void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, - const T &m, - X &theta, - bool &unlimited) { - const X &x = this->m_x[j]; - const X &lbound = this->m_lower_bounds[j]; - - if (this->below_bound(x, lbound)) { - limit_theta((lbound - x) / m, theta, unlimited); - } else { - const X &ubound = this->m_upper_bounds[j]; - if (this->below_bound(x, ubound)) { - limit_theta((ubound - x) / m, theta, unlimited); - } else if (!this->above_bound(x, ubound)) { - theta = zero_of_type(); - unlimited = false; - } - } - } - - void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, - const T &m, - X &theta, - bool &unlimited) { - // lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed); - const X &x = this->m_x[j]; - const X &ubound = this->m_upper_bounds[j]; - if (this->above_bound(x, ubound)) { - limit_theta((ubound - x) / m, theta, unlimited); - } else { - const X &lbound = this->m_lower_bounds[j]; - if (this->above_bound(x, lbound)) { - limit_theta((lbound - x) / m, theta, unlimited); - } else if (!this->below_bound(x, lbound)) { - theta = zero_of_type(); - unlimited = false; - } - } - } - - void limit_theta_on_basis_column_for_feas_case_m_pos_no_check( - unsigned j, const T &m, X &theta, bool &unlimited) { - lp_assert(m > 0); - limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited); - if (theta < zero_of_type()) { - theta = zero_of_type(); - } - } - - // j is a basic column or the entering, in any case x[j] has to stay feasible. - // m is the multiplier. updating t in a way that holds the following - // x[j] + t * m >= this->m_lower_bounds[j]( if m < 0 ) - // or - // x[j] + t * m <= this->m_upper_bounds[j] ( if m > 0) - void limit_theta_on_basis_column(unsigned j, T m, X &theta, bool &unlimited) { - switch (this->m_column_types[j]) { - case column_type::free_column: - break; - case column_type::upper_bound: - if (this->current_x_is_feasible()) { - if (m > 0) - limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, - unlimited); - } else { // inside of feasibility_loop - if (m > 0) - limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( - j, m, theta, unlimited); - else - limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( - j, m, theta, unlimited); - } - break; - case column_type::lower_bound: - if (this->current_x_is_feasible()) { - if (m < 0) - limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, - unlimited); - } else { - if (m < 0) - limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( - j, m, theta, unlimited); - else - limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( - j, m, theta, unlimited); - } - break; - // case fixed: - // if (get_this->current_x_is_feasible()) { - // theta = zero_of_type(); - // break; - // } - // if (m < 0) - // limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m, - // theta); - // else - // limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m, - // theta); - // break; - case column_type::fixed: - case column_type::boxed: - if (this->current_x_is_feasible()) { - if (m > 0) { - limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, - unlimited); - } else { - limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, - unlimited); + return true; } - } else { - if (m > 0) { - limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta, - unlimited); - } else { - limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta, - unlimited); + + void limit_inf_on_lower_bound_m_pos(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets larger + lp_assert(m > 0); + if (this->below_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } } - } - break; - default: - UNREACHABLE(); - } - if (!unlimited && theta < zero_of_type()) { - theta = zero_of_type(); - } - } + void limit_inf_on_upper_bound_m_neg(const T &m, const X &x, const X &bound, + X &theta, bool &unlimited) { + // x gets smaller + lp_assert(m < 0); + if (this->above_bound(x, bound)) { + limit_theta((bound - x) / m, theta, unlimited); + } + } - bool column_is_benefitial_for_entering_basis(unsigned j) const; - void init_infeasibility_costs(); - void print_column(unsigned j, std::ostream &out); + void limit_theta_on_basis_column_for_inf_case_m_pos_boxed(unsigned j, + const T &m, + X &theta, + bool &unlimited) { + const X &x = this->m_x[j]; + const X &lbound = this->m_lower_bounds[j]; - void print_bound_info_and_x(unsigned j, std::ostream &out); + if (this->below_bound(x, lbound)) { + limit_theta((lbound - x) / m, theta, unlimited); + } else { + const X &ubound = this->m_upper_bounds[j]; + if (this->below_bound(x, ubound)) { + limit_theta((ubound - x) / m, theta, unlimited); + } else if (!this->above_bound(x, ubound)) { + theta = zero_of_type(); + unlimited = false; + } + } + } - bool basis_column_is_set_correctly(unsigned j) const { - return this->m_A.m_columns[j].size() == 1; - } + void limit_theta_on_basis_column_for_inf_case_m_neg_boxed(unsigned j, + const T &m, + X &theta, + bool &unlimited) { + // lp_assert(m < 0 && this->m_column_type[j] == column_type::boxed); + const X &x = this->m_x[j]; + const X &ubound = this->m_upper_bounds[j]; + if (this->above_bound(x, ubound)) { + limit_theta((ubound - x) / m, theta, unlimited); + } else { + const X &lbound = this->m_lower_bounds[j]; + if (this->above_bound(x, lbound)) { + limit_theta((lbound - x) / m, theta, unlimited); + } else if (!this->below_bound(x, lbound)) { + theta = zero_of_type(); + unlimited = false; + } + } + } - bool basis_columns_are_set_correctly() const { - for (unsigned j : this->m_basis) - if (!basis_column_is_set_correctly(j)) - return false; + void limit_theta_on_basis_column_for_feas_case_m_pos_no_check( + unsigned j, const T &m, X &theta, bool &unlimited) { + lp_assert(m > 0); + limit_theta((this->m_upper_bounds[j] - this->m_x[j]) / m, theta, unlimited); + if (theta < zero_of_type()) { + theta = zero_of_type(); + } + } - return this->m_basis_heading.size() == this->m_A.column_count() && - this->m_basis.size() == this->m_A.row_count(); - } + // j is a basic column or the entering, in any case x[j] has to stay feasible. + // m is the multiplier. updating t in a way that holds the following + // x[j] + t * m >= this->m_lower_bounds[j]( if m < 0 ) + // or + // x[j] + t * m <= this->m_upper_bounds[j] ( if m > 0) + void limit_theta_on_basis_column(unsigned j, T m, X &theta, bool &unlimited) { + switch (this->m_column_types[j]) { + case column_type::free_column: + break; + case column_type::upper_bound: + if (this->current_x_is_feasible()) { + if (m > 0) + limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, + unlimited); + } else { // inside of feasibility_loop + if (m > 0) + limit_theta_on_basis_column_for_inf_case_m_pos_upper_bound( + j, m, theta, unlimited); + else + limit_theta_on_basis_column_for_inf_case_m_neg_upper_bound( + j, m, theta, unlimited); + } + break; + case column_type::lower_bound: + if (this->current_x_is_feasible()) { + if (m < 0) + limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, + unlimited); + } else { + if (m < 0) + limit_theta_on_basis_column_for_inf_case_m_neg_lower_bound( + j, m, theta, unlimited); + else + limit_theta_on_basis_column_for_inf_case_m_pos_lower_bound( + j, m, theta, unlimited); + } + break; + // case fixed: + // if (get_this->current_x_is_feasible()) { + // theta = zero_of_type(); + // break; + // } + // if (m < 0) + // limit_theta_on_basis_column_for_inf_case_m_neg_fixed(j, m, + // theta); + // else + // limit_theta_on_basis_column_for_inf_case_m_pos_fixed(j, m, + // theta); + // break; + case column_type::fixed: + case column_type::boxed: + if (this->current_x_is_feasible()) { + if (m > 0) { + limit_theta_on_basis_column_for_feas_case_m_pos_no_check(j, m, theta, + unlimited); + } else { + limit_theta_on_basis_column_for_feas_case_m_neg_no_check(j, m, theta, + unlimited); + } + } else { + if (m > 0) { + limit_theta_on_basis_column_for_inf_case_m_pos_boxed(j, m, theta, + unlimited); + } else { + limit_theta_on_basis_column_for_inf_case_m_neg_boxed(j, m, theta, + unlimited); + } + } - void init_run_tableau(); - void update_x_tableau(unsigned entering, const X &delta); - // the delta is between the old and the new cost (old - new) - void update_reduced_cost_for_basic_column_cost_change(const T &delta, - unsigned j) { - lp_assert(this->m_basis_heading[j] >= 0); - unsigned i = static_cast(this->m_basis_heading[j]); - for (const row_cell &rc : this->m_A.m_rows[i]) { - unsigned k = rc.var(); - if (k == j) - continue; - this->m_d[k] += delta * rc.coeff(); - } - } + break; + default: + UNREACHABLE(); + } + if (!unlimited && theta < zero_of_type()) { + theta = zero_of_type(); + } + } - bool update_basis_and_x_tableau(int entering, int leaving, X const &tt); - void init_reduced_costs_tableau(); - void init_tableau_rows() { - m_bland_mode_tableau = false; - m_left_basis_tableau.clear(); - m_left_basis_tableau.resize(this->m_A.column_count()); - m_left_basis_repeated = 0; - } - // stage1 constructor - lp_primal_core_solver( - static_matrix &A, - 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, const vector &column_type_array, - const vector &lower_bound_values, const 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_bland_mode_threshold(1000) { - this->set_status(lp_status::UNKNOWN); - } + bool column_is_benefitial_for_entering_basis(unsigned j) const; + void init_infeasibility_costs(); + void print_column(unsigned j, std::ostream &out); - friend core_solver_pretty_printer; -}; + void print_bound_info_and_x(unsigned j, std::ostream &out); + + bool basis_column_is_set_correctly(unsigned j) const { + return this->m_A.m_columns[j].size() == 1; + } + + bool basis_columns_are_set_correctly() const { + for (unsigned j : this->m_basis) + if (!basis_column_is_set_correctly(j)) + return false; + + return this->m_basis_heading.size() == this->m_A.column_count() && + this->m_basis.size() == this->m_A.row_count(); + } + + void init_run_tableau(); + void update_x_tableau(unsigned entering, const X &delta); + // the delta is between the old and the new cost (old - new) + void update_reduced_cost_for_basic_column_cost_change(const T &delta, + unsigned j) { + lp_assert(this->m_basis_heading[j] >= 0); + unsigned i = static_cast(this->m_basis_heading[j]); + for (const row_cell &rc : this->m_A.m_rows[i]) { + unsigned k = rc.var(); + if (k == j) + continue; + this->m_d[k] += delta * rc.coeff(); + } + } + + bool update_basis_and_x_tableau(int entering, int leaving, X const &tt); + void init_reduced_costs_tableau(); + void init_tableau_rows() { + m_bland_mode_tableau = false; + m_left_basis_tableau.clear(); + m_left_basis_tableau.resize(this->m_A.column_count()); + m_left_basis_repeated = 0; + } + // stage1 constructor + lp_primal_core_solver( + static_matrix &A, + 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, const vector &column_type_array, + const vector &lower_bound_values, const 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_bland_mode_threshold(1000) { + this->set_status(lp_status::UNKNOWN); + } + + friend core_solver_pretty_printer; + }; } // namespace lp diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h index 7ae8624af..68b062d64 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -162,13 +162,14 @@ template void lp_primal_core_solver::advance_on_en return; } if (!is_zero(t)) { - if (this->current_x_is_feasible() ) { + if (this->current_x_is_feasible()) { if (m_sign_of_entering_delta == -1) t = -t; } this->update_basis_and_x_tableau(entering, leaving, t); this->iters_with_no_cost_growing() = 0; - } else { + } + else { this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); this->change_basis(entering, leaving); } diff --git a/src/smt/theory_arith_int.h b/src/smt/theory_arith_int.h index 496d27614..911cd2eca 100644 --- a/src/smt/theory_arith_int.h +++ b/src/smt/theory_arith_int.h @@ -908,15 +908,20 @@ namespace smt { inf_numeral l, u; numeral m; unsigned num_patched = 0, num_not_patched = 0; + unsigned num_one = 0, num_divides = 0; for (theory_var v = 0; v < num; v++) { if (!is_non_base(v)) continue; get_freedom_interval(v, inf_l, l, inf_u, u, m); - if (m.is_one() && get_value(v).is_int()) + if (m.is_one() && get_value(v).is_int()) { + num_one++; continue; + } // check whether value of v is already a multiple of m. - if ((get_value(v).get_rational() / m).is_int()) + if ((get_value(v).get_rational() / m).is_int()) { + num_divides++; continue; + } TRACE("patch_int", tout << "TARGET v" << v << " -> ["; @@ -927,13 +932,14 @@ namespace smt { tout << ", m: " << m << ", val: " << get_value(v) << ", is_int: " << is_int(v) << "\n";); +#if 0 verbose_stream() << "TARGET v" << v << " -> ["; if (inf_l) verbose_stream() << "-oo"; else verbose_stream() << ceil(l); verbose_stream() << ", "; if (inf_u) verbose_stream() << "oo"; else verbose_stream() << floor(u); verbose_stream() << "]"; verbose_stream() << ", m: " << m << ", val: " << get_value(v) << ", is_int: " << is_int(v) << "\n"; - +#endif if (!inf_l) l = ceil(l); if (!inf_u) @@ -957,7 +963,9 @@ namespace smt { set_value(v, inf_numeral(0)); } SASSERT(m_to_patch.empty()); - verbose_stream() << "patched " << num_patched << " not patched " << num_not_patched << "\n"; + verbose_stream() << "patched " << num_patched << " not patched " << num_not_patched << " ones " << num_one << " divides " << num_divides << "\n"; + //display(verbose_stream()); + //exit(0); } /**