From 16b71fe911a706c3b14d495bca3577ff58a5c247 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Wed, 18 Jul 2018 12:50:54 -0700 Subject: [PATCH 1/7] work on static_matrix's cells Signed-off-by: Lev trying the new scheme in static_matrix : in progress Signed-off-by: Lev Nachmanson in the middle of changes in static_matrix Signed-off-by: Lev Nachmanson more fixes in static_matrix.h Signed-off-by: Lev Nachmanson debug Signed-off-by: Lev Nachmanson fixes in static_matrix Signed-off-by: Lev fixes in static_matrix, column_strip Signed-off-by: Lev fixes in static_matrix Signed-off-by: Lev Nachmanson fixes for static_matrix Signed-off-by: Lev work on static_matrix Signed-off-by: Lev work on static_matrix Signed-off-by: Lev progress in static_matrix Signed-off-by: Lev fix a bug in swap_with_head_cell Signed-off-by: Lev progress in static_matrix Signed-off-by: Lev compress rows and columns if needed Signed-off-by: Lev fix in compression of cells Signed-off-by: Lev Nachmanson --- src/shell/main.cpp | 8 +- src/tactic/smtlogics/qflia_tactic.cpp | 4 + src/test/lp/lp.cpp | 2 +- src/util/lp/bound_analyzer_on_row.h | 70 +-- src/util/lp/bound_propagator.cpp | 2 + src/util/lp/core_solver_pretty_printer_def.h | 2 +- src/util/lp/indexed_vector.h | 1 + src/util/lp/int_solver.cpp | 17 +- src/util/lp/lar_core_solver.h | 10 +- src/util/lp/lar_core_solver_def.h | 5 +- src/util/lp/lar_solver.cpp | 66 ++- src/util/lp/lp_core_solver_base.h | 18 +- src/util/lp/lp_core_solver_base_def.h | 91 +-- src/util/lp/lp_primal_core_solver.h | 62 +- src/util/lp/lp_primal_core_solver_def.h | 2 +- .../lp/lp_primal_core_solver_tableau_def.h | 22 +- src/util/lp/lu.h | 4 +- src/util/lp/lu_def.h | 2 +- src/util/lp/random_updater_def.h | 5 +- src/util/lp/static_matrix.cpp | 9 +- src/util/lp/static_matrix.h | 549 +++++++++++++++--- src/util/lp/static_matrix_def.h | 204 +++---- 22 files changed, 767 insertions(+), 388 deletions(-) diff --git a/src/shell/main.cpp b/src/shell/main.cpp index d367b8ec6..f6c809bb2 100644 --- a/src/shell/main.cpp +++ b/src/shell/main.cpp @@ -295,7 +295,13 @@ void parse_cmd_line_args(int argc, char ** argv) { int STD_CALL main(int argc, char ** argv) { - try{ +#ifdef DUMP_ARGS + std::cout << "args are: "; + for (int i = 0; i < argc; i++) + std::cout << argv[i] <<" "; + std::cout << std::endl; +#endif + try{ unsigned return_value = 0; memory::initialize(0); memory::exit_when_out_of_memory(true, "ERROR: out of memory"); diff --git a/src/tactic/smtlogics/qflia_tactic.cpp b/src/tactic/smtlogics/qflia_tactic.cpp index eed4e4425..ad5b660c3 100644 --- a/src/tactic/smtlogics/qflia_tactic.cpp +++ b/src/tactic/smtlogics/qflia_tactic.cpp @@ -212,6 +212,9 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) { tactic * st = using_params(and_then(preamble_st, +#if 1 + mk_smt_tactic()), +#else or_else(mk_ilp_model_finder_tactic(m), mk_pb_tactic(m), and_then(fail_if_not(mk_is_quasi_pb_probe()), @@ -219,6 +222,7 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) { mk_fail_if_undecided_tactic()), mk_bounded_tactic(m), mk_smt_tactic())), +#endif main_p); diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 6e418fe68..37a096827 100644 --- a/src/test/lp/lp.cpp +++ b/src/test/lp/lp.cpp @@ -2541,7 +2541,7 @@ void read_row_cols(unsigned i, static_matrix& A, std::ifstream & lp_assert(r.size() == 4); unsigned j = atoi(r[1].c_str()); double v = atof(r[3].c_str()); - A.set(i, j, v); + A.add_new_element(i, j, v); } while (true); } diff --git a/src/util/lp/bound_analyzer_on_row.h b/src/util/lp/bound_analyzer_on_row.h index 196551f20..3d047cd86 100644 --- a/src/util/lp/bound_analyzer_on_row.h +++ b/src/util/lp/bound_analyzer_on_row.h @@ -29,65 +29,7 @@ Revision History: namespace lp { template // C plays a role of a container class bound_analyzer_on_row { - struct term_with_basis_col { - const C & m_row; - unsigned m_bj; - struct ival { - unsigned m_var; - const mpq & m_coeff; - ival(unsigned var, const mpq & val) : m_var(var), m_coeff(val) { - } - unsigned var() const { return m_var;} - const mpq & coeff() const { return m_coeff; } - }; - - term_with_basis_col(const C& row, unsigned bj) : m_row(row), m_bj(bj) {} - struct const_iterator { - // fields - typename C::const_iterator m_it; - unsigned m_bj; - - - //typedefs - - - typedef const_iterator self_type; - typedef ival value_type; - typedef ival reference; - typedef int difference_type; - typedef std::forward_iterator_tag iterator_category; - - reference operator*() const { - if (m_bj == static_cast(-1)) - return ival((*m_it).var(), (*m_it).coeff()); - return ival(m_bj, - 1); - } - self_type operator++() { self_type i = *this; operator++(1); return i; } - - self_type operator++(int) { - if (m_bj == static_cast(-1)) - m_it++; - else - m_bj = static_cast(-1); - return *this; - } - - // constructor - const_iterator(const typename C::const_iterator& it, unsigned bj) : - m_it(it), - m_bj(bj) - {} - bool operator==(const self_type &other) const { - return m_it == other.m_it && m_bj == other.m_bj ; - } - bool operator!=(const self_type &other) const { return !(*this == other); } - }; - const_iterator begin() const { - return const_iterator( m_row.begin(), m_bj); - } - const_iterator end() const { return const_iterator(m_row.end(), m_bj); } - }; - term_with_basis_col m_row; + const C& m_row; bound_propagator & m_bp; unsigned m_row_or_term_index; int m_column_of_u; // index of an unlimited from above monoid @@ -105,7 +47,7 @@ public : bound_propagator & bp ) : - m_row(it, bj), + m_row(it), m_bp(bp), m_row_or_term_index(row_or_term_index), m_column_of_u(-1), @@ -117,6 +59,8 @@ public : unsigned j; void analyze() { for (const auto & c : m_row) { + if (c.dead()) + continue; if ((m_column_of_l == -2) && (m_column_of_u == -2)) break; analyze_bound_on_var_on_coeff(c.var(), c.coeff()); @@ -226,6 +170,7 @@ public : mpq total; lp_assert(is_zero(total)); for (const auto& p : m_row) { + if (p.dead()) continue; bool str; total -= monoid_min(p.coeff(), p.var(), str); if (str) @@ -234,6 +179,7 @@ public : for (const auto &p : m_row) { + if (p.dead()) continue; bool str; bool a_is_pos = is_pos(p.coeff()); mpq bound = total / p.coeff() + monoid_min_no_mult(a_is_pos, p.var(), str); @@ -251,6 +197,7 @@ public : mpq total; lp_assert(is_zero(total)); for (const auto &p : m_row) { + if (p.dead()) continue; bool str; total -= monoid_max(p.coeff(), p.var(), str); if (str) @@ -258,6 +205,7 @@ public : } for (const auto& p : m_row) { + if (p.dead()) continue; bool str; bool a_is_pos = is_pos(p.coeff()); mpq bound = total / p.coeff() + monoid_max_no_mult(a_is_pos, p.var(), str); @@ -280,6 +228,7 @@ public : mpq bound = -m_rs.x; bool strict = false; for (const auto& p : m_row) { + if (p.dead()) continue; j = p.var(); if (j == static_cast(m_column_of_u)) { u_coeff = p.coeff(); @@ -309,6 +258,7 @@ public : mpq bound = -m_rs.x; bool strict = false; for (const auto &p : m_row) { + if (p.dead()) continue; j = p.var(); if (j == static_cast(m_column_of_l)) { l_coeff = p.coeff(); diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index c4fa2aefa..199fef9b0 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -16,6 +16,8 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const { return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]; } void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) { + TRACE("try_add_bound", + tout << "v = " << v << ", j = " << j << std::endl;); j = m_lar_solver.adjust_column_index_to_term_index(j); if (m_lar_solver.is_term(j)) { // lp treats terms as not having a free coefficient, restoring it below for the outside consumption diff --git a/src/util/lp/core_solver_pretty_printer_def.h b/src/util/lp/core_solver_pretty_printer_def.h index 58ffbb481..dcf3e0a11 100644 --- a/src/util/lp/core_solver_pretty_printer_def.h +++ b/src/util/lp/core_solver_pretty_printer_def.h @@ -100,7 +100,7 @@ template void core_solver_pretty_printer::init_m_ for (unsigned column = 0; column < ncols(); column++) { vector t(nrows(), zero_of_type()); for (const auto & c : m_core_solver.m_A.m_columns[column]){ - t[c.m_i] = m_core_solver.m_A.get_val(c); + t[c.var()] = m_core_solver.m_A.get_val(c); } string name = m_core_solver.column_name(column); diff --git a/src/util/lp/indexed_vector.h b/src/util/lp/indexed_vector.h index d472a4a1e..24078ce84 100644 --- a/src/util/lp/indexed_vector.h +++ b/src/util/lp/indexed_vector.h @@ -172,6 +172,7 @@ public: } unsigned var() const { return m_var;} const T & coeff() const { return m_coeff; } + bool dead() const { return false; } }; struct const_iterator { diff --git a/src/util/lp/int_solver.cpp b/src/util/lp/int_solver.cpp index 2e744e3d6..78206dc53 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -100,6 +100,7 @@ bool int_solver::is_gomory_cut_target(const row_strip& row) { // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). unsigned j; for (const auto & p : row) { + if (p.dead()) continue; j = p.var(); if (is_base(j)) continue; if (!at_bound(j)) @@ -311,6 +312,7 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row TRACE("gomory_cut", tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl; for (auto & p : row) { + if (p.dead()) continue; m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout); } tout << "inf_col = " << inf_col << std::endl; @@ -324,7 +326,8 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row bool some_int_columns = false; mpq f_0 = int_solver::fractional_part(get_value(inf_col)); mpq one_min_f_0 = 1 - f_0; - for (auto & p : row) { + for (const auto & p : row) { + if (p.dead()) continue; x_j = p.var(); if (x_j == inf_col) continue; @@ -353,6 +356,7 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row int int_solver::find_free_var_in_gomory_row(const row_strip& row) { unsigned j; for (const auto & p : row) { + if (p.dead()) continue; j = p.var(); if (!is_base(j) && is_free(j)) return static_cast(j); @@ -789,6 +793,7 @@ lia_move int_solver::patch_nbasic_columns() { mpq get_denominators_lcm(const row_strip & row) { mpq r(1); for (auto & c : row) { + if (c.dead()) continue; r = lcm(r, denominator(c.coeff())); } return r; @@ -802,6 +807,7 @@ bool int_solver::gcd_test_for_row(static_matrix> & A, uns bool least_coeff_is_bounded = false; unsigned j; for (auto &c : A.m_rows[i]) { + if (c.dead()) continue; j = c.var(); const mpq& a = c.coeff(); if (m_lar_solver->column_is_fixed(j)) { @@ -867,6 +873,7 @@ void int_solver::add_to_explanation_from_fixed_or_boxed_column(unsigned j) { } void int_solver::fill_explanation_from_fixed_columns(const row_strip & row) { for (const auto & c : row) { + if (c.dead()) continue; if (!m_lar_solver->column_is_fixed(c.var())) continue; add_to_explanation_from_fixed_or_boxed_column(c.var()); @@ -892,6 +899,7 @@ bool int_solver::ext_gcd_test(const row_strip & row, mpq a; unsigned j; for (const auto & c : row) { + if (c.dead()) continue; j = c.var(); const mpq & a = c.coeff(); if (m_lar_solver->column_is_fixed(j)) @@ -1023,6 +1031,7 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq lp_assert(settings().use_tableau()); const auto & A = m_lar_solver->A_r(); for (const auto &c : A.column(j)) { + if (c.dead()) continue; row_index = c.var(); const mpq & a = c.coeff(); @@ -1152,13 +1161,15 @@ bool int_solver::at_upper(unsigned j) const { void int_solver::display_row_info(std::ostream & out, unsigned row_index) const { auto & rslv = m_lar_solver->m_mpq_lar_core_solver.m_r_solver; - for (auto &c: rslv.m_A.m_rows[row_index]) { + for (const auto &c: rslv.m_A.m_rows[row_index]) { + if (c.dead()) continue; if (numeric_traits::is_pos(c.coeff())) out << "+"; out << c.coeff() << rslv.column_name(c.var()) << " "; } - for (auto& c: rslv.m_A.m_rows[row_index]) { + for (const auto& c: rslv.m_A.m_rows[row_index]) { + if (c.dead()) continue; rslv.print_column_bound_info(c.var(), out); } rslv.print_column_bound_info(rslv.m_basis[row_index], out); diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h index 5f5518528..4223f5964 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -339,7 +339,7 @@ public: if (!update_xj_and_get_delta(j, pos_type, delta)) continue; for (const auto & cc : m_r_solver.m_A.m_columns[j]){ - unsigned i = cc.m_i; + unsigned i = cc.var(); unsigned jb = m_r_solver.m_basis[i]; m_r_solver.update_x_with_delta_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc)); } @@ -583,8 +583,7 @@ public: if (!m_r_solver.m_settings.use_tableau()) return true; for (unsigned j : m_r_solver.m_basis) { - lp_assert(m_r_solver.m_A.m_columns[j].size() == 1); - lp_assert(m_r_solver.m_A.get_val(m_r_solver.m_A.m_columns[j][0]) == one_of_type()); + lp_assert(m_r_solver.m_A.m_columns[j].live_size() == 1); } for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) { if (m_r_solver.m_basis_heading[j] >= 0) continue; @@ -633,8 +632,9 @@ public: void create_double_matrix(static_matrix & A) { for (unsigned i = 0; i < m_r_A.row_count(); i++) { auto & row = m_r_A.m_rows[i]; - for (row_cell & c : row) { - A.set(i, c.m_j, c.get_val().get_double()); + for (row_cell & c : row.m_cells) { + if (c.dead()) continue; + A.add_new_element(i, c.var(), c.get_val().get_double()); } } } diff --git a/src/util/lp/lar_core_solver_def.h b/src/util/lp/lar_core_solver_def.h index f9937e77a..c626362fc 100644 --- a/src/util/lp/lar_core_solver_def.h +++ b/src/util/lp/lar_core_solver_def.h @@ -226,8 +226,9 @@ void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() { unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau]; m_infeasible_sum_sign = m_r_solver.inf_sign_of_column(bj); m_infeasible_linear_combination.clear(); - for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) { - m_infeasible_linear_combination.push_back(std::make_pair( rc.get_val(), rc.m_j)); + for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau].m_cells) { + if (rc.dead()) continue; + m_infeasible_linear_combination.push_back(std::make_pair( rc.get_val(), rc.var())); } } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index c4f7ff196..5408bd922 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -162,9 +162,8 @@ void lar_solver::analyze_new_bounds_on_row_tableau( unsigned row_index, bound_propagator & bp ) { - if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation) + if (A_r().m_rows[row_index].live_size() > settings().max_row_length_for_bound_propagation) return; - lp_assert(use_tableau()); bound_analyzer_on_row>::analyze_row(A_r().m_rows[row_index], static_cast(-1), @@ -216,7 +215,8 @@ void lar_solver::explain_implied_bound(implied_bound & ib, bound_propagator & bp bound_j = m_var_register.external_to_local(bound_j); } for (auto const& r : A_r().m_rows[i]) { - unsigned j = r.m_j; + if (r.dead()) continue; + unsigned j = r.var(); if (j == bound_j) continue; mpq const& a = r.get_val(); int a_sign = is_pos(a)? 1: -1; @@ -428,8 +428,9 @@ void lar_solver::set_costs_to_zero(const lar_term& term) { if (i < 0) jset.insert(j); else { - for (auto & rc : A_r().m_rows[i]) - jset.insert(rc.m_j); + for (const auto & rc : A_r().m_rows[i]) + if (rc.alive()) + jset.insert(rc.var()); } } @@ -639,7 +640,8 @@ void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column(unsigned j void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j) { for (auto & rc : m_mpq_lar_core_solver.m_r_A.m_columns[j]) - m_rows_with_changed_bounds.insert(rc.m_i); + if (rc.alive()) + m_rows_with_changed_bounds.insert(rc.var()); } bool lar_solver::use_tableau() const { return m_settings.use_tableau(); } @@ -667,8 +669,10 @@ void lar_solver::adjust_x_of_column(unsigned j) { bool lar_solver::row_is_correct(unsigned i) const { numeric_pair r = zero_of_type>(); - for (const auto & c : A_r().m_rows[i]) - r += c.m_value * m_mpq_lar_core_solver.m_r_x[c.m_j]; + for (const auto & c : A_r().m_rows[i]) { + if (c.dead()) continue; + r += c.coeff() * m_mpq_lar_core_solver.m_r_x[c.var()]; + } return is_zero(r); } @@ -691,7 +695,8 @@ bool lar_solver::costs_are_used() const { void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair & delta) { if (use_tableau()) { for (const auto & c : A_r().m_columns[j]) { - unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.m_i]; + if (c.dead()) continue; + unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()]; if (tableau_with_costs()) { m_basic_columns_with_changed_cost.insert(bj); } @@ -790,10 +795,10 @@ numeric_pair lar_solver::get_basic_var_value_from_row_directly(unsigned i) unsigned bj = m_mpq_lar_core_solver.m_r_solver.m_basis[i]; for (const auto & c: A_r().m_rows[i]) { - if (c.m_j == bj) continue; - const auto & x = m_mpq_lar_core_solver.m_r_x[c.m_j]; + if (c.var() == bj) continue; + const auto & x = m_mpq_lar_core_solver.m_r_x[c.var()]; if (!is_zero(x)) - r -= c.m_value * x; + r -= c.coeff() * x; } return r; } @@ -866,14 +871,14 @@ void lar_solver::fill_last_row_of_A_r(static_matrix> & A, lp_assert(A.row_count() > 0); lp_assert(A.column_count() > 0); unsigned last_row = A.row_count() - 1; - lp_assert(A.m_rows[last_row].size() == 0); + lp_assert(A.m_rows[last_row].live_size() == 0); for (auto & t : ls->m_coeffs) { lp_assert(!is_zero(t.second)); var_index j = t.first; - A.set(last_row, j, - t.second); + A.add_new_element(last_row, j, - t.second); } unsigned basis_j = A.column_count() - 1; - A.set(last_row, basis_j, mpq(1)); + A.add_new_element(last_row, basis_j, mpq(1)); } template @@ -895,7 +900,7 @@ void lar_solver::copy_from_mpq_matrix(static_matrix & matr) { matr.m_columns.resize(A_r().column_count()); for (unsigned i = 0; i < matr.row_count(); i++) { for (auto & it : A_r().m_rows[i]) { - matr.set(i, it.m_j, convert_struct::convert(it.get_val())); + matr.add_new_element(i, it.var(), convert_struct::convert(it.get_val())); } } } @@ -1327,16 +1332,17 @@ void lar_solver::make_sure_that_the_bottom_right_elem_not_zero_in_tableau(unsign lp_assert(A_r().row_count() == i + 1 && A_r().column_count() == j + 1); auto & last_column = A_r().m_columns[j]; int non_zero_column_cell_index = -1; - for (unsigned k = last_column.size(); k-- > 0;){ + for (unsigned k = last_column.cells_size(); k-- > 0;){ auto & cc = last_column[k]; - if (cc.m_i == i) + if (cc.dead()) continue; + if (cc.var() == i) return; non_zero_column_cell_index = k; } lp_assert(non_zero_column_cell_index != -1); lp_assert(static_cast(non_zero_column_cell_index) != i); - m_mpq_lar_core_solver.m_r_solver.transpose_rows_tableau(last_column[non_zero_column_cell_index].m_i, i); + m_mpq_lar_core_solver.m_r_solver.transpose_rows_tableau(last_column[non_zero_column_cell_index].var(), i); } void lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { @@ -1351,24 +1357,29 @@ void lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { auto & last_row = A_r().m_rows[i]; mpq &cost_j = m_mpq_lar_core_solver.m_r_solver.m_costs[j]; bool cost_is_nz = !is_zero(cost_j); - for (unsigned k = last_row.size(); k-- > 0;) { + for (unsigned k = last_row.cells_size(); k-- > 0;) { auto &rc = last_row[k]; + if (rc.dead()) { + last_row.pop(); + continue; + } if (cost_is_nz) { - m_mpq_lar_core_solver.m_r_solver.m_d[rc.m_j] += cost_j*rc.get_val(); + m_mpq_lar_core_solver.m_r_solver.m_d[rc.var()] += cost_j*rc.get_val(); } - A_r().remove_element(last_row, rc); + A_r().remove_element(rc); } - lp_assert(last_row.size() == 0); - lp_assert(A_r().m_columns[j].size() == 0); + lp_assert(last_row.live_size() == 0); + lp_assert(A_r().m_columns[j].live_size() == 0); A_r().m_rows.pop_back(); A_r().m_columns.pop_back(); + lp_assert(A_r().is_correct()); slv.m_b.pop_back(); } void lar_solver::remove_last_column_from_A() { // the last column has to be empty - lp_assert(A_r().m_columns.back().size() == 0); + lp_assert(A_r().m_columns.back().live_size() == 0); A_r().m_columns.pop_back(); } @@ -1840,11 +1851,12 @@ void lar_solver::fill_last_row_of_A_d(static_matrix & A, const l for (auto & t : ls->m_coeffs) { lp_assert(!is_zero(t.second)); var_index j = t.first; - A.set(last_row, j, -t.second.get_double()); + A.add_new_element(last_row, j, -t.second.get_double()); } unsigned basis_j = A.column_count() - 1; - A.set(last_row, basis_j, -1); + A.add_new_element(last_row, basis_j, -1); + lp_assert(A.is_correct()); } void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) { diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index b8b7b7c0d..35a989c07 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -213,8 +213,17 @@ public: unsigned m = m_A.row_count(); for (unsigned i = 0; i < m; i++) { unsigned bj = m_basis[i]; - lp_assert(m_A.m_columns[bj].size() > 0); - if (m_A.m_columns[bj].size() > 1 || m_A.get_val(m_A.m_columns[bj][0]) != one_of_type()) return true; + lp_assert(m_A.m_columns[bj].live_size() > 0); + if (m_A.m_columns[bj].live_size() > 1) + return true; + for (const auto & c : m_A.m_columns[bj]) { + if (c.dead()) + continue; + if (m_A.get_val(c) != one_of_type()) + return true; + else + break; + } } return false; } @@ -237,8 +246,9 @@ public: } } else { auto d = m_costs[j]; - for (auto & cc : this->m_A.m_columns[j]) { - d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc); + for (const auto & cc : this->m_A.m_columns[j]) { + if (cc.dead()) continue; + d -= this->m_costs[this->m_basis[cc.var()]] * this->m_A.get_val(cc); } if (m_d[j] != d) { return false; diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index f20eaf042..ddae7bf20 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -103,9 +103,12 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) { T &a = m_d[j]; if (is_zero(a)) return; - for (const row_cell & r: m_A.m_rows[i]) - if (r.m_j != j) - m_d[r.m_j] -= a * r.get_val(); + for (const row_cell & r: m_A.m_rows[i]){ + if (r.dead()) continue; + if (r.var() != j) + m_d[r.var()] -= a * r.get_val(); + } +ls a = zero_of_type(); // zero the pivot column's m_d finally } @@ -308,8 +311,9 @@ calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned pivot_row) { if (numeric_traits::is_zero(pi_1)) { continue; } - for (auto & c : m_A.m_rows[i]) { - unsigned j = c.m_j; + for (auto & c : m_A.m_rows[i].m_cells) { + if (c.dead()) continue; + unsigned j = c.var(); if (m_basis_heading[j] < 0) { m_pivot_row.add_value_at_index_with_drop_tolerance(j, c.get_val() * pi_1); } @@ -331,7 +335,7 @@ update_x(unsigned entering, const X& delta) { } else for (const auto & c : m_A.m_columns[entering]) { - unsigned i = c.m_i; + unsigned i = c.var(); m_x[m_basis[i]] -= delta * m_A.get_val(c); } } @@ -449,7 +453,7 @@ rs_minus_Anx(vector & rs) { while (row--) { auto &rsv = rs[row] = m_b[row]; for (auto & it : m_A.m_rows[row]) { - unsigned j = it.m_j; + unsigned j = it.var(); if (m_basis_heading[j] < 0) { rsv -= m_x[j] * it.get_val(); } @@ -589,9 +593,11 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) { lp_assert(numeric_traits::precise()); int pivot_index = -1; auto & row = m_A.m_rows[pivot_row]; - unsigned size = row.size(); + unsigned size = row.cells_size(); for (unsigned j = 0; j < size; j++) { - if (row[j].m_j == pivot_col) { + auto & c = row[j]; + if (c.dead()) continue; + if (c.var() == pivot_col) { pivot_index = static_cast(j); break; } @@ -599,26 +605,33 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) { if (pivot_index == -1) return false; auto & pivot_cell = row[pivot_index]; - if (is_zero(pivot_cell.m_value)) + T & coeff = pivot_cell.coeff(); + if (is_zero(coeff)) return false; - this->m_b[pivot_row] /= pivot_cell.m_value; + this->m_b[pivot_row] /= coeff; for (unsigned j = 0; j < size; j++) { - if (row[j].m_j != pivot_col) { - row[j].m_value /= pivot_cell.m_value; + auto & c = row[j]; + if (c.dead()) continue; + if (c.var() != pivot_col) { + c.coeff() /= coeff; } } - pivot_cell.m_value = one_of_type(); + coeff = one_of_type(); + lp_assert(m_A.is_correct()); return true; } template bool lp_core_solver_base:: pivot_column_tableau(unsigned j, unsigned piv_row_index) { - if (!divide_row_by_pivot(piv_row_index, j)) + lp_assert(m_A.is_correct()); + m_A.compress_row_if_needed(piv_row_index); + if (!divide_row_by_pivot(piv_row_index, j)) return false; auto &column = m_A.m_columns[j]; int pivot_col_cell_index = -1; - for (unsigned k = 0; k < column.size(); k++) { - if (column[k].m_i == piv_row_index) { + for (unsigned k = 0; k < column.cells_size(); k++) { + if (column[k].dead()) continue; + if (column[k].var() == piv_row_index) { pivot_col_cell_index = k; break; } @@ -626,25 +639,26 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { if (pivot_col_cell_index < 0) return false; - if (pivot_col_cell_index != 0) { - lp_assert(column.size() > 1); - // swap the pivot column cell with the head cell - auto c = column[0]; - column[0] = column[pivot_col_cell_index]; - column[pivot_col_cell_index] = c; + if (pivot_col_cell_index != 0) + m_A.swap_with_head_cell(j, pivot_col_cell_index); - m_A.m_rows[piv_row_index][column[0].m_offset].m_offset = 0; - m_A.m_rows[c.m_i][c.m_offset].m_offset = pivot_col_cell_index; - } - while (column.size() > 1) { + lp_assert(m_A.is_correct()); + while (column.live_size() > 1) { auto & c = column.back(); - lp_assert(c.m_i != piv_row_index); + if (c.dead()) { + column.pop_last_dead_cell(); + continue; + } + lp_assert(c.var() != piv_row_index); if(! m_A.pivot_row_to_row_given_cell(piv_row_index, c, j)) { return false; } if (m_pivoted_rows!= nullptr) - m_pivoted_rows->insert(c.m_i); + m_pivoted_rows->insert(c.var()); } + m_A.compress_column_if_needed(j); + lp_assert(column.live_size() == 1); + lp_assert(m_A.is_correct()); if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs) pivot_to_reduced_costs_tableau(piv_row_index, j); @@ -758,10 +772,11 @@ fill_reduced_costs_from_m_y_by_rows() { while (i--) { const T & y = m_y[i]; if (is_zero(y)) continue; - for (row_cell & it : m_A.m_rows[i]) { - j = it.m_j; + for (row_cell & c : m_A.m_rows[i].m_cells) { + if (c.dead()) continue; + j = c.var(); if (m_basis_heading[j] < 0) { - m_d[j] -= y * it.get_val(); + m_d[j] -= y * c.get_val(); } } } @@ -802,7 +817,7 @@ find_error_in_BxB(vector& rs){ while (row--) { auto &rsv = rs[row]; for (auto & it : m_A.m_rows[row]) { - unsigned j = it.m_j; + unsigned j = it.var(); if (m_basis_heading[j] >= 0) { rsv -= m_x[j] * it.get_val(); } @@ -957,7 +972,8 @@ template void lp_core_solver_base::pivot_fixed_v if (get_column_type(basic_j) != column_type::fixed) continue; T a; unsigned j; - for (auto &c : m_A.m_rows[i]) { + for (auto &c : m_A.m_rows[i].m_cells) { + if (c.dead()) continue; j = c.var(); if (j == basic_j) continue; @@ -972,7 +988,8 @@ template void lp_core_solver_base::pivot_fixed_v template bool lp_core_solver_base::remove_from_basis(unsigned basic_j) { indexed_vector w(m_basis.size()); // the buffer unsigned i = m_basis_heading[basic_j]; - for (auto &c : m_A.m_rows[i]) { + for (auto &c : m_A.m_rows[i].m_cells) { + if (c.dead()) continue; if (c.var() == basic_j) continue; if (pivot_column_general(c.var(), basic_j, w)) @@ -1043,8 +1060,8 @@ void lp_core_solver_base::calculate_pivot_row(unsigned i) { if (m_settings.use_tableau()) { unsigned basic_j = m_basis[i]; for (auto & c : m_A.m_rows[i]) { - if (c.m_j != basic_j) - m_pivot_row.set_value(c.get_val(), c.m_j); + if (c.var() != basic_j) + m_pivot_row.set_value(c.get_val(), c.var()); } return; } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 26fc550b3..2700fc168 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -172,7 +172,7 @@ public: bool monoid_can_decrease(const row_cell & rc) const { - unsigned j = rc.m_j; + unsigned j = rc.var(); lp_assert(this->column_is_feasible(j)); switch (this->m_column_types[j]) { case column_type::free_column: @@ -205,7 +205,7 @@ public: } bool monoid_can_increase(const row_cell & rc) const { - unsigned j = rc.m_j; + unsigned j = rc.var(); lp_assert(this->column_is_feasible(j)); switch (this->m_column_types[j]) { case column_type::free_column: @@ -239,8 +239,9 @@ public: unsigned get_number_of_basic_vars_that_might_become_inf(unsigned j) const { // consider looking at the signs here: todo unsigned r = 0; - for (auto & cc : this->m_A.m_columns[j]) { - unsigned k = this->m_basis[cc.m_i]; + for (const auto & cc : this->m_A.m_columns[j]) { + if (cc.dead()) continue; + unsigned k = this->m_basis[cc.var()]; if (this->m_column_types[k] != column_type::free_column) r++; } @@ -253,7 +254,8 @@ public: 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.m_j == bj) + if (rc.dead()) continue; + if (rc.var() == bj) continue; if (bj_needs_to_grow) { if (!monoid_can_decrease(rc)) @@ -262,9 +264,9 @@ public: if (!monoid_can_increase(rc)) continue; } - if (rc.m_j < static_cast(j) ) { - j = rc.m_j; - a_ent = rc.m_value; + if (rc.var() < static_cast(j) ) { + j = rc.var(); + a_ent = rc.coeff(); } } if (j == -1) { @@ -278,13 +280,18 @@ public: 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 - vector*> choices; - unsigned num_of_non_free_basics = 1000000; + int choice = -1; + int nchoices = 0; unsigned len = 100000000; 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]) { - unsigned j = rc.m_j; + for (unsigned k = 0; k < this->m_A.m_rows[i].m_cells.size(); k++) { + const row_cell& rc = this->m_A.m_rows[i].m_cells[k]; + if (rc.dead()) continue; + const row_cell& rc = this->m_A.m_rows[i].m_cells[k]; + if (rc.dead()) continue; +>>>>>>> e6c612f... trying the new scheme in static_matrix : in progress + unsigned j = rc.var(); if (j == bj) continue; if (bj_needs_to_grow) { @@ -297,26 +304,24 @@ public: 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(); - choices.clear(); - choices.push_back(&rc); + len = this->m_A.m_columns[j].live_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() % 2)) { - choices.push_back(&rc); - len = this->m_A.m_columns[j].size(); + this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % (++nchoices))) { + choice = k; + len = this->m_A.m_columns[j].live_size(); } } - if (choices.size() == 0) { + if (choice == -1) { m_inf_row_index_for_tableau = i; return -1; } - const row_cell* rc = choices.size() == 1? choices[0] : - choices[this->m_settings.random_next() % choices.size()]; - - a_ent = rc->m_value; - return rc->m_j; + const row_cell& rc = this->m_A.m_rows[i].m_cells[choice]; + a_ent = rc.coeff(); + return rc.var(); } static X positive_infinity() { return convert_struct::convert(std::numeric_limits::max()); @@ -827,11 +832,11 @@ public: this->m_rows_nz.resize(this->m_A.row_count()); for (unsigned i = 0; i < this->m_A.column_count(); i++) { if (this->m_columns_nz[i] == 0) - this->m_columns_nz[i] = this->m_A.m_columns[i].size(); + this->m_columns_nz[i] = this->m_A.m_columns[i].live_size(); } for (unsigned i = 0; i < this->m_A.row_count(); i++) { if (this->m_rows_nz[i] == 0) - this->m_rows_nz[i] = this->m_A.m_rows[i].size(); + this->m_rows_nz[i] = this->m_A.m_rows[i].live_size(); } } @@ -861,7 +866,7 @@ public: unsigned solve_with_tableau(); bool basis_column_is_set_correctly(unsigned j) const { - return this->m_A.m_columns[j].size() == 1; + return this->m_A.m_columns[j].live_size() == 1; } @@ -881,7 +886,8 @@ public: 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.m_j; + if (rc.dead()) continue; + unsigned k = rc.var(); if (k == j) continue; this->m_d[k] += delta * rc.get_val(); diff --git a/src/util/lp/lp_primal_core_solver_def.h b/src/util/lp/lp_primal_core_solver_def.h index 1e9edbd31..eb4a9ce1d 100644 --- a/src/util/lp/lp_primal_core_solver_def.h +++ b/src/util/lp/lp_primal_core_solver_def.h @@ -975,7 +975,7 @@ template void lp_primal_core_solver::delete_fa template void lp_primal_core_solver::init_column_norms() { lp_assert(numeric_traits::precise() == false); for (unsigned j = 0; j < this->m_n(); j++) { - this->m_column_norms[j] = T(static_cast(this->m_A.m_columns[j].size() + 1)) + this->m_column_norms[j] = T(static_cast(this->m_A.m_columns[j].live_size() + 1)) + T(static_cast(this->m_settings.random_next() % 10000)) / T(100000); } diff --git a/src/util/lp/lp_primal_core_solver_tableau_def.h b/src/util/lp/lp_primal_core_solver_tableau_def.h index f85c131a3..be28bf36f 100644 --- a/src/util/lp/lp_primal_core_solver_tableau_def.h +++ b/src/util/lp/lp_primal_core_solver_tableau_def.h @@ -266,17 +266,18 @@ template int lp_primal_core_solver::find_leaving_ unsigned row_min_nz = this->m_n() + 1; m_leaving_candidates.clear(); auto & col = this->m_A.m_columns[entering]; - unsigned col_size = col.size(); + unsigned col_size = col.cells_size(); for (;k < col_size && unlimited; k++) { const column_cell & c = col[k]; - unsigned i = c.m_i; + if (c.dead()) continue; + unsigned i = c.var(); const T & ed = this->m_A.get_val(c); lp_assert(!numeric_traits::is_zero(ed)); unsigned j = this->m_basis[i]; limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, t, unlimited); if (!unlimited) { m_leaving_candidates.push_back(j); - row_min_nz = this->m_A.m_rows[i].size(); + row_min_nz = this->m_A.m_rows[i].live_size(); } } if (unlimited) { @@ -288,14 +289,15 @@ template int lp_primal_core_solver::find_leaving_ X ratio; for (;k < col_size; k++) { const column_cell & c = col[k]; - unsigned i = c.m_i; + if (c.dead()) continue; + unsigned i = c.var(); const T & ed = this->m_A.get_val(c); lp_assert(!numeric_traits::is_zero(ed)); unsigned j = this->m_basis[i]; unlimited = true; limit_theta_on_basis_column(j, -ed * m_sign_of_entering_delta, ratio, unlimited); if (unlimited) continue; - unsigned i_nz = this->m_A.m_rows[i].size(); + unsigned i_nz = this->m_A.m_rows[i].live_size(); if (ratio < t) { t = ratio; m_leaving_candidates.clear(); @@ -304,7 +306,7 @@ template int lp_primal_core_solver::find_leaving_ } else if (ratio == t && i_nz < row_min_nz) { m_leaving_candidates.clear(); m_leaving_candidates.push_back(j); - row_min_nz = this->m_A.m_rows[i].size(); + row_min_nz = this->m_A.m_rows[i].live_size(); } else if (ratio == t && i_nz == row_min_nz) { m_leaving_candidates.push_back(j); } @@ -348,6 +350,7 @@ template void lp_primal_core_solver::init_run_tab template bool lp_primal_core_solver:: update_basis_and_x_tableau(int entering, int leaving, X const & tt) { lp_assert(this->use_tableau()); + lp_assert(entering != leaving); update_x_tableau(entering, tt); this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); this->change_basis(entering, leaving); @@ -358,7 +361,8 @@ update_x_tableau(unsigned entering, const X& delta) { this->add_delta_to_x_and_call_tracker(entering, delta); if (!this->m_using_infeas_costs) { for (const auto & c : this->m_A.m_columns[entering]) { - unsigned i = c.m_i; + if (c.dead()) continue; + unsigned i = c.var(); this->update_x_with_delta_and_track_feasibility(this->m_basis[i], - delta * this->m_A.get_val(c)); } } else { // m_using_infeas_costs == true @@ -366,7 +370,7 @@ update_x_tableau(unsigned entering, const X& delta) { lp_assert(this->m_costs[entering] == zero_of_type()); // m_d[entering] can change because of the cost change for basic columns. for (const auto & c : this->m_A.m_columns[entering]) { - unsigned i = c.m_i; + unsigned i = c.var(); unsigned j = this->m_basis[i]; this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c)); update_inf_cost_for_column_tableau(j); @@ -407,7 +411,7 @@ template void lp_primal_core_solver::init_reduced else { T& d = this->m_d[j] = this->m_costs[j]; for (auto & cc : this->m_A.m_columns[j]) { - d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc); + d -= this->m_costs[this->m_basis[cc.var()]] * this->m_A.get_val(cc); } } } diff --git a/src/util/lp/lu.h b/src/util/lp/lu.h index 3b5a9cc59..95c00414e 100644 --- a/src/util/lp/lu.h +++ b/src/util/lp/lu.h @@ -330,11 +330,11 @@ public: for (unsigned i = m_prev; i < m; i++) { for (const row_cell & c : m_A.m_rows[i]) { - int h = heading[c.m_j]; + int h = heading[c.var()]; if (h < 0) { continue; } - columns_to_replace.insert(c.m_j); + columns_to_replace.insert(c.var()); } } return columns_to_replace; diff --git a/src/util/lp/lu_def.h b/src/util/lp/lu_def.h index 51f291e7e..00d0466b7 100644 --- a/src/util/lp/lu_def.h +++ b/src/util/lp/lu_def.h @@ -386,7 +386,7 @@ void lu< M>::find_error_of_yB_indexed(const indexed_vector& y, const vector::add_columns_at_the_end(unsigned int); +template void static_matrix::add_new_element(unsigned i, unsigned j, const double & v); +template void static_matrix::add_new_element(unsigned i, unsigned j, const mpq & v); +template void static_matrix::add_new_element(unsigned i, unsigned j, const mpq & v); template void static_matrix::clear(); #ifdef Z3DEBUG template bool static_matrix::is_correct() const; @@ -47,7 +50,6 @@ template double static_matrix::get_min_abs_in_row(unsigned int) template void static_matrix::init_empty_matrix(unsigned int, unsigned int); template void static_matrix::init_row_columns(unsigned int, unsigned int); template static_matrix::ref & static_matrix::ref::operator=(double const&); -template void static_matrix::set(unsigned int, unsigned int, double const&); template static_matrix::static_matrix(unsigned int, unsigned int); template void static_matrix::add_column_to_vector(mpq const&, unsigned int, mpq*) const; template void static_matrix::add_columns_at_the_end(unsigned int); @@ -63,7 +65,6 @@ template mpq static_matrix::get_min_abs_in_column(unsigned int) const; template mpq static_matrix::get_min_abs_in_row(unsigned int) const; template void static_matrix::init_row_columns(unsigned int, unsigned int); template static_matrix::ref& static_matrix::ref::operator=(mpq const&); -template void static_matrix::set(unsigned int, unsigned int, mpq const&); template static_matrix::static_matrix(unsigned int, unsigned int); #ifdef Z3DEBUG @@ -72,13 +73,11 @@ template bool static_matrix >::is_correct() const; template void static_matrix >::copy_column_to_indexed_vector(unsigned int, indexed_vector&) const; template mpq static_matrix >::get_elem(unsigned int, unsigned int) const; template void static_matrix >::init_empty_matrix(unsigned int, unsigned int); -template void static_matrix >::set(unsigned int, unsigned int, mpq const&); - template bool lp::static_matrix::pivot_row_to_row_given_cell(unsigned int, column_cell &, unsigned int); template bool lp::static_matrix::pivot_row_to_row_given_cell(unsigned int, column_cell& , unsigned int); template bool lp::static_matrix >::pivot_row_to_row_given_cell(unsigned int, column_cell&, unsigned int); -template void lp::static_matrix >::remove_element(vector, true, unsigned int>&, lp::row_cell&); +template void lp::static_matrix >::remove_element(lp::row_cell&); } diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index e4bf257f4..036f628b7 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -1,22 +1,22 @@ /*++ -Copyright (c) 2017 Microsoft Corporation + Copyright (c) 2017 Microsoft Corporation -Module Name: + Module Name: - + -Abstract: + Abstract: - + -Author: + Author: - Lev Nachmanson (levnach) + Lev Nachmanson (levnach) -Revision History: + Revision History: ---*/ + --*/ #pragma once #include "util/vector.h" @@ -29,30 +29,387 @@ Revision History: #include namespace lp { -struct column_cell { - unsigned m_i; // points to the row - unsigned m_offset; // the offset of the element in the matrix row - column_cell(unsigned i, unsigned offset) : m_i(i), m_offset(offset) { +class column_cell { + unsigned m_i; // points to the row + unsigned m_offset; // the offset of the element in the matrix row, or the next dead cell in the column_strip +public: + column_cell(unsigned i, unsigned offset) : m_i(i), m_offset(offset) { } + column_cell(unsigned i) : m_i(i) { } + + // index of of the cell row + unsigned var() const { + lp_assert(alive()); + return m_i; + } + unsigned &var() { + return m_i; + } + unsigned offset() const { + lp_assert(alive()); + return m_offset; + } + unsigned & offset() { + lp_assert(alive()); + return m_offset; + } + + unsigned next_dead_index() const { + lp_assert(dead()); + return m_offset; + } + unsigned & next_dead_index() { + return m_offset; } + bool alive() const { return !dead(); } + bool dead() const { return m_i == static_cast(-1); } + void set_dead() { m_i = -1;} + + + }; template -struct row_cell { - unsigned m_j; // points to the column - unsigned m_offset; // offset in column - T m_value; +class row_cell { + unsigned m_j; // points to the column + unsigned m_offset; // offset in column, or the next dead cell in the row_strip + T m_value; +public: row_cell(unsigned j, unsigned offset, T const & val) : m_j(j), m_offset(offset), m_value(val) { } - const T & get_val() const { return m_value;} - T & get_val() { return m_value;} - void set_val(const T& v) { m_value = v; } - unsigned var() const { return m_j;} - const T & coeff() const { return m_value;} + row_cell(unsigned j, T const & val) : m_j(j), m_value(val) { + } + const T & get_val() const { + lp_assert(alive()); + return m_value; + } + T & get_val() { + lp_assert(alive()); + return m_value; + } + void set_val(const T& v) { + lp_assert(alive()); + m_value = v; + } + // index of the cell column + unsigned var() const { + lp_assert(alive()); + return m_j; + } + unsigned &var() { + return m_j; + } + const T & coeff() const { + lp_assert(alive()); + return m_value; + } + T & coeff() { + lp_assert(alive()); + return m_value; + } + // offset in the column vector + unsigned offset() const { + lp_assert(alive()); + return m_offset; + } + unsigned & offset() { + return m_offset; + } + unsigned next_dead_index() const { + lp_assert(dead()); + return m_offset; + } + unsigned & next_dead_index() { + lp_assert(dead()); + return m_offset; + } + bool alive() const { return !dead(); } + bool dead() const { return m_j == static_cast(-1); } + void set_dead() { m_j = static_cast(-1); } + + }; template -using row_strip = vector>; +class row_strip { + unsigned m_live_size; + int m_first_dead; +public: + row_strip() : m_live_size(0), m_first_dead(-1) {} + row_cell* begin() { return m_cells.begin();} + const row_cell* begin() const { return m_cells.begin();} + row_cell* end() { return m_cells.end();} + const row_cell* end() const { return m_cells.end();} + unsigned live_size() const { return m_live_size; } + vector> m_cells; + unsigned cells_size() const { return m_cells.size(); } + const row_cell & operator[](unsigned i) const { return m_cells[i]; } + row_cell & operator[](unsigned i) { return m_cells[i];} + + void skip_dead_cell(unsigned k) { + lp_assert(m_cells[k].dead()) + auto * c = &m_cells[m_first_dead]; + for (; c->next_dead_index() != k; c = &m_cells[c->next_dead_index()]); + lp_assert(c->next_dead_index() == k); + c->next_dead_index() = m_cells[k].next_dead_index(); + } + + void pop_last_dead_cell() { + lp_assert(m_cells.back().dead()); + unsigned last = m_cells.size() - 1; + if (m_first_dead == static_cast(last)) + m_first_dead = m_cells.back().next_dead_index(); + else { + skip_dead_cell(last); + } + m_cells.pop_back(); + } + + void pop(){ + bool dead = m_cells.back().dead(); + if (dead) { + pop_last_dead_cell(); + } else { + m_live_size --; + m_cells.pop_back(); + } + } + + bool empty() const { return m_live_size == 0; } + + void delete_at(unsigned j) { + lp_assert(m_cells[j].alive()); + m_live_size--; + if (j == m_cells.size() - 1) + m_cells.pop_back(); + else { + auto & c = m_cells[j]; + c.set_dead(); + c.next_dead_index() = m_first_dead; + m_first_dead = j; + } + lp_assert(is_correct()); + } + + bool is_correct() const { + std::set d0; + std::set d1; + unsigned alive = 0; + for (unsigned i = 0; i < m_cells.size(); i++) { + if (m_cells[i].dead()) + d0.insert(i); + else + alive ++; + } + + if ( alive != m_live_size) + return false; + + for (unsigned i = m_first_dead; i != static_cast(-1) ; i = m_cells[i].next_dead_index()) + d1.insert(i); + return d0 == d1; + } + + row_cell & add_cell(unsigned j, const T& val, unsigned & index) { +#ifdef Z3DEBUG + for (const auto & c : m_cells) { + if (c.dead()) continue; + if (c.var() == j) + std::cout << "oops" << std::endl; + } +#endif + if (m_first_dead != -1) { + auto & ret = m_cells[index = m_first_dead]; + m_first_dead = ret.next_dead_index(); + m_live_size++; + ret.var() = j; + ret.coeff() = val; + return ret; + } + lp_assert(m_live_size == m_cells.size()); + index = m_live_size++; + m_cells.push_back(row_cell(j, val)); + return m_cells.back(); + } + + void shrink_to_live() { + m_cells.shrink(live_size()); + m_first_dead = -1; + } +}; + +class column_strip { + vector m_cells; + unsigned m_live_size; + int m_first_dead; +public: + column_strip() : m_live_size(0), m_first_dead(-1) { } + column_cell* begin() { return m_cells.begin();} + const column_cell* begin() const { return m_cells.begin();} + column_cell* end() { return m_cells.end();} + const column_cell* end() const { return m_cells.end();} + unsigned live_size() const { + return m_live_size; + } + + unsigned cells_size() const { + return m_cells.size(); + } + + unsigned cells_size() const { + return m_cells.size(); + } + + column_cell& back() { return m_cells.back(); } + void delete_at(unsigned j) { + lp_assert(m_cells[j].alive()); + m_live_size--; + if (j == m_cells.size() - 1) + m_cells.pop_back(); + else { + auto & c = m_cells[j]; + c.set_dead(); + c.next_dead_index() = m_first_dead; + m_first_dead = j; + } + lp_assert(is_correct()); + } + + const column_cell& operator[] (unsigned i) const { return m_cells[i];} + column_cell& operator[](unsigned i) { return m_cells[i];} + void pop() { + bool dead = m_cells.back().dead(); + if (dead) { + pop_last_dead_cell(); + } else { + m_live_size --; + m_cells.pop_back(); + } + } + void skip_dead_cell(unsigned k) { + lp_assert(m_cells[k].dead()) + auto * c = &m_cells[m_first_dead]; + for (; c->next_dead_index() != k; c = &m_cells[c->next_dead_index()]); + lp_assert(c->next_dead_index() == k); + c->next_dead_index() = m_cells[k].next_dead_index(); + } + + void pop_last_dead_cell() { + lp_assert(m_cells.back().dead()); + unsigned last = m_cells.size() - 1; + if (m_first_dead == static_cast(last)) + m_first_dead = m_cells.back().next_dead_index(); + else { + skip_dead_cell(last); + } + m_cells.pop_back(); + } + + std::set d0; + std::set d1; + unsigned alive = 0; + for (unsigned i = 0; i < m_cells.size(); i++) { + if (m_cells[i].dead()) + d0.insert(i); + else + alive ++; + } + + if (alive != m_live_size) + return false; + for (unsigned i = m_first_dead; i != static_cast(-1) ; i = m_cells[i].next_dead_index()) + d1.insert(i); + return d0 == d1; + } + + void swap_with_head_cell(unsigned i) { + lp_assert(i > 0); + lp_assert(m_cells[i].alive()); + column_cell head_copy = m_cells[0]; + if (head_copy.dead()) { + if (m_first_dead == 0) { + m_first_dead = i; + } else { + column_cell * c = &m_cells[m_first_dead]; + for (; c->next_dead_index() != 0; c = &m_cells[c->next_dead_index()]); + lp_assert(c->next_dead_index() == 0); + c->next_dead_index() = i; + } + } + m_cells[0] = m_cells[i]; + m_cells[i] = head_copy; + lp_assert(is_correct()); + } + + void swap_with_head_cell(unsigned i) { + lp_assert(i > 0); + lp_assert(m_cells[i].alive()); + column_cell head_copy = m_cells[0]; + if (head_copy.dead()) { + if (m_first_dead == 0) { + m_first_dead = i; + } else { + column_cell * c = &m_cells[m_first_dead]; + for (; c->next_dead_index() != 0; c = &m_cells[c->next_dead_index()]); + lp_assert(c->next_dead_index() == 0); + c->next_dead_index() = i; + } + } + m_cells[0] = m_cells[i]; + m_cells[i] = head_copy; + lp_assert(is_correct()); + } + + column_cell & add_cell(unsigned i, unsigned & index) { + if (m_first_dead != -1) { + auto & ret = m_cells[index = m_first_dead]; + m_first_dead = ret.next_dead_index(); + m_live_size++; + ret.var() = i; + return ret; + } + lp_assert(m_live_size == m_cells.size()); + index = m_live_size++; + m_cells.push_back(column_cell(i)); + return m_cells.back(); + } + column_cell & add_cell(unsigned i, unsigned & index) { + if (m_first_dead != -1) { + auto & ret = m_cells[index = m_first_dead]; + m_first_dead = ret.next_dead_index(); + m_live_size++; + ret.var() = i; + return ret; + } + lp_assert(m_live_size == m_cells.size()); + index = m_live_size++; + m_cells.push_back(column_cell(i)); + return m_cells.back(); + } + void shrink_to_live() { + m_cells.shrink(live_size()); + m_first_dead = -1; + } +}; + +template +void compress_cells(A& cells, vector& vec_of_cells) { + if (2 * cells.live_size() < cells.cells_size()) + return; + unsigned j = 0; // the available target for copy + for (unsigned i = 0; i < cells.cells_size(); i++) { + auto & c = cells[i]; + if (c.dead()) continue; + if (i == j) { + j++; + continue; + } + vec_of_cells[c.var()][c.offset()].offset() = j; + cells[j++] = c; + } + cells.shrink_to_live(); +} + // each assignment for this matrix should be issued only once!!! template @@ -66,13 +423,13 @@ class static_matrix unsigned m_n; dim(unsigned m, unsigned n) :m_m(m), m_n(n) {} }; - std::stack m_stack; + // fields + std::stack m_stack; public: - typedef vector column_strip; - vector m_vector_of_row_offsets; - indexed_vector m_work_vector; - vector> m_rows; - vector m_columns; + vector m_vector_of_row_offsets; + indexed_vector m_work_vector; + vector> m_rows; + vector m_columns; // starting inner classes class ref { static_matrix & m_matrix; @@ -80,9 +437,9 @@ public: unsigned m_col; public: ref(static_matrix & m, unsigned row, unsigned col):m_matrix(m), m_row(row), m_col(col) {} - ref & operator=(T const & v) { m_matrix.set( m_row, m_col, v); return *this; } + ref & operator=(T const & v) { m_matrix.add_new_element( m_row, m_col, v); return *this; } - ref operator=(ref & v) { m_matrix.set(m_row, m_col, v.m_matrix.get(v.m_row, v.m_col)); return *this; } + ref operator=(ref & v) { m_matrix.add_new_element(m_row, m_col, v.m_matrix.get(v.m_row, v.m_col)); return *this; } operator T () const { return m_matrix.get_elem(m_row, m_col); } }; @@ -98,7 +455,7 @@ public: public: const T & get_val(const column_cell & c) const { - return m_rows[c.m_i][c.m_offset].get_val(); + return m_rows[c.var()][c.offset()].get_val(); } column_cell & get_column_cell(const row_cell &rc) { @@ -107,7 +464,7 @@ public: void init_row_columns(unsigned m, unsigned n); - // constructor with no parameters + // constructor with no parameters static_matrix() {} // constructor @@ -142,12 +499,12 @@ public: void remove_last_column(unsigned j); - void remove_element(vector> & row, row_cell & elem_to_remove); + void remove_element(row_cell & elem_to_remove); void multiply_column(unsigned column, T const & alpha) { for (auto & t : m_columns[column]) { - auto & r = m_rows[t.m_i][t.m_offset]; - r.m_value *= alpha; + auto & r = m_rows[t.var()].m_cells[t.offset()]; + r.coeff() *= alpha; } } @@ -165,8 +522,6 @@ public: return column_cell(column, offset); } - void set(unsigned row, unsigned col, T const & val); - ref operator()(unsigned row, unsigned col) { return ref(*this, row, col); } std::set> get_domain(); @@ -175,8 +530,8 @@ public: T get_max_abs_in_row(unsigned row) const; void add_column_to_vector (const T & a, unsigned j, T * v) const { - for (const auto & it : m_columns[j]) { - v[it.m_i] += a * get_val(it); + for (const auto & c : m_columns[j]) { + v[c.var()] += a * get_val(c); } } @@ -189,9 +544,6 @@ public: void check_consistency(); #endif - - void cross_out_row(unsigned k); - // void fix_row_indices_in_each_column_for_crossed_row(unsigned k); @@ -202,9 +554,9 @@ public: T get_elem(unsigned i, unsigned j) const; - unsigned number_of_non_zeroes_in_column(unsigned j) const { return m_columns[j].size(); } + unsigned number_of_non_zeroes_in_column(unsigned j) const { return m_columns[j].live_size(); } - unsigned number_of_non_zeroes_in_row(unsigned i) const { return m_rows[i].size(); } + unsigned number_of_non_zeroes_in_row(unsigned i) const { return m_rows[i].live_size(); } unsigned number_of_non_zeroes() const { unsigned ret = 0; @@ -237,12 +589,12 @@ public: m_stack.push(d); } - void pop_row_columns(const vector> & row) { + void pop_row_columns(const row_strip & row) { for (auto & c : row) { - unsigned j = c.m_j; - auto & col = m_columns[j]; - lp_assert(col[col.size() - 1].m_i == m_rows.size() -1 ); // todo : start here!!!! - col.pop_back(); + if (c.dead()) { + continue; + } + m_columns[c.var()].delete_at(c.offset()); } } @@ -272,8 +624,9 @@ public: } void multiply_row(unsigned row, T const & alpha) { - for (auto & t : m_rows[row]) { - t.m_value *= alpha; + for (auto & t : m_rows[row].m_cells) { + if (t.alive()) + t.coeff() *= alpha; } } @@ -286,8 +639,8 @@ public: T dot_product_with_column(const vector & y, unsigned j) const { lp_assert(j < column_count()); T ret = numeric_traits::zero(); - for (auto & it : m_columns[j]) { - ret += y[it.m_i] * get_val(it); // get_value_of_column_cell(it); + for (auto & c : m_columns[j]) { + ret += y[c.var()] * get_val(c); } return ret; } @@ -301,18 +654,20 @@ public: m_rows[i] = m_rows[ii]; m_rows[ii] = t; // now fix the columns - for (auto & rc : m_rows[i]) { - column_cell & cc = m_columns[rc.m_j][rc.m_offset]; - lp_assert(cc.m_i == ii); - cc.m_i = i; + for (const auto & rc : m_rows[i]) { + if (rc.dead()) continue; + column_cell & cc = m_columns[rc.var()][rc.offset()]; + lp_assert(cc.var() == ii); + cc.var() = i; } - for (auto & rc : m_rows[ii]) { - column_cell & cc = m_columns[rc.m_j][rc.m_offset]; - lp_assert(cc.m_i == i); - cc.m_i = ii; + for (const auto & rc : m_rows[ii]) { + if (rc.dead()) continue; + column_cell & cc = m_columns[rc.var()][rc.offset()]; + lp_assert(cc.var() == i); + cc.var() = ii; } - } + void fill_last_row_with_pivoting_loop_block(unsigned j, const vector & basis_heading) { int row_index = basis_heading[j]; if (row_index < 0) @@ -322,17 +677,18 @@ public: return; for (const auto & c : m_rows[row_index]) { - if (c.m_j == j) { + if (c.dead()) continue; + if (c.var() == j) { continue; } - T & wv = m_work_vector.m_data[c.m_j]; + T & wv = m_work_vector.m_data[c.var()]; bool was_zero = is_zero(wv); - wv -= alpha * c.m_value; + wv -= alpha * c.coeff(); if (was_zero) - m_work_vector.m_index.push_back(c.m_j); + m_work_vector.m_index.push_back(c.var()); else { if (is_zero(wv)) { - m_work_vector.erase_from_index(c.m_j); + m_work_vector.erase_from_index(c.var()); } } } @@ -350,7 +706,7 @@ public: lp_assert(row_count() > 0); m_work_vector.resize(column_count()); T a; - // we use the form -it + 1 = 0 + // we use the form -it + 1 = 0 m_work_vector.set_value(one_of_type(), bj); for (auto p : row) { m_work_vector.set_value(-p.coeff(), p.var()); @@ -366,10 +722,10 @@ public: unsigned last_row = row_count() - 1; for (unsigned j : m_work_vector.m_index) { - set (last_row, j, m_work_vector.m_data[j]); + add_new_element(last_row, j, m_work_vector.m_data[j]); } lp_assert(column_count() > 0); - set(last_row, column_count() - 1, one_of_type()); + add_new_element(last_row, column_count() - 1, one_of_type()); } void copy_column_to_vector (unsigned j, vector & v) const { @@ -377,7 +733,7 @@ public: for (auto & it : m_columns[j]) { const T& val = get_val(it); if (!is_zero(val)) - v[it.m_i] = val; + v[it.var()] = val; } } @@ -386,7 +742,8 @@ public: L ret = zero_of_type(); lp_assert(row < m_rows.size()); for (auto & it : m_rows[row]) { - ret += w[it.m_j] * it.get_val(); + if (it.dead()) continue; + ret += w[it.var()] * it.coeff(); } return ret; } @@ -397,9 +754,9 @@ public: // constructor column_cell_plus(const column_cell & c, const static_matrix& A) : m_c(c), m_A(A) {} - unsigned var() const { return m_c.m_i; } - const T & coeff() const { return m_A.m_rows[var()][m_c.m_offset].get_val(); } - + unsigned var() const { return m_c.var(); } + const T & coeff() const { return m_A.m_rows[var()][m_c.offset()].get_val(); } + bool dead() const { return m_c.dead(); } }; struct column_container { @@ -411,7 +768,7 @@ public: // fields const column_cell *m_c; const static_matrix& m_A; - + const column_cell *m_end; //typedefs @@ -425,13 +782,19 @@ public: reference operator*() const { return column_cell_plus(*m_c, m_A); } - self_type operator++() { self_type i = *this; m_c++; return i; } - self_type operator++(int) { m_c++; return *this; } + self_type operator++() { self_type i = *this; + m_c++; + return i; + } + + self_type operator++(int) { + m_c++; + return *this; + } const_iterator(const column_cell* it, const static_matrix& A) : m_c(it), - m_A(A) - {} + m_A(A){} bool operator==(const self_type &other) const { return m_c == other.m_c; } @@ -449,8 +812,32 @@ public: column_container column(unsigned j) const { return column_container(j, *this); + } + void swap_with_head_cell(unsigned j, unsigned offset) { + column_strip & col = m_columns[j]; + column_cell & head = col[0]; + column_cell & cc = col[offset]; + if (head.alive()) { + m_rows[head.var()][head.offset()].offset() = offset; + } + lp_assert(cc.alive()); + m_rows[cc.var()][cc.offset()].offset() = 0; + col.swap_with_head_cell(offset); + } + + + void compress_row_if_needed(unsigned i) { + compress_cells(m_rows[i], m_columns); + lp_assert(is_correct()); + } + + void compress_column_if_needed(unsigned j) { + compress_cells(m_columns[j], m_rows); + lp_assert(is_correct()); + } + ref_row operator[](unsigned i) const { return ref_row(*this, i);} typedef T coefftype; typedef X argtype; diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index 42249b4d8..698363289 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -36,44 +36,54 @@ void static_matrix::init_row_columns(unsigned m, unsigned n) { template void static_matrix::scan_row_ii_to_offset_vector(const row_strip & rvals) { - for (unsigned j = 0; j < rvals.size(); j++) - m_vector_of_row_offsets[rvals[j].m_j] = j; + for (unsigned j = 0; j < rvals.cells_size(); j++) { + if (rvals[j].dead()) continue; + m_vector_of_row_offsets[rvals[j].var()] = j; + } } template bool static_matrix::pivot_row_to_row_given_cell(unsigned i, column_cell & c, unsigned pivot_col) { - unsigned ii = c.m_i; + lp_assert(is_correct()); + unsigned ii = c.var(); lp_assert(i < row_count() && ii < column_count() && i != ii); T alpha = -get_val(c); lp_assert(!is_zero(alpha)); + compress_row_if_needed(ii); auto & rowii = m_rows[ii]; - remove_element(rowii, rowii[c.m_offset]); + remove_element(rowii[c.offset()]); scan_row_ii_to_offset_vector(rowii); - unsigned prev_size_ii = rowii.size(); + unsigned prev_size_ii = rowii.cells_size(); // run over the pivot row and update row ii - for (const auto & iv : m_rows[i]) { - unsigned j = iv.m_j; + for (const auto & iv : m_rows[i].m_cells) { + if (iv.dead()) continue; + unsigned j = iv.var(); if (j == pivot_col) continue; - T alv = alpha * iv.m_value; - lp_assert(!is_zero(iv.m_value)); + lp_assert(!is_zero(iv.coeff())); + T alv = alpha * iv.coeff(); int j_offs = m_vector_of_row_offsets[j]; if (j_offs == -1) { // it is a new element add_new_element(ii, j, alv); } else { - rowii[j_offs].m_value += alv; + rowii[j_offs].coeff() += alv; } } // clean the work vector for (unsigned k = 0; k < prev_size_ii; k++) { - m_vector_of_row_offsets[rowii[k].m_j] = -1; + auto & c = rowii[k]; + if (c.dead()) continue; + m_vector_of_row_offsets[c.var()] = -1; } - // remove zeroes - for (unsigned k = rowii.size(); k-- > 0; ) { - if (is_zero(rowii[k].m_value)) - remove_element(rowii, rowii[k]); + for (unsigned k = rowii.cells_size(); k-- > 0; ) { + auto & c = rowii[k]; + if (c.dead()) + continue; + if (is_zero(c.coeff())) + remove_element(c); } + lp_assert(is_correct()); return !rowii.empty(); } @@ -86,7 +96,7 @@ static_matrix::static_matrix(static_matrix const &A, unsigned * /* basis * init_row_columns(m, m); while (m--) { for (auto & col : A.m_columns[m]){ - set(col.m_i, m, A.get_value_of_column_cell(col)); + set(col.var(), m, A.get_value_of_column_cell(col)); } } } @@ -107,14 +117,14 @@ template void static_matrix::init_empty_matrix init_row_columns(m, n); } -template unsigned static_matrix::lowest_row_in_column(unsigned col) { - lp_assert(col < column_count()); +template unsigned static_matrix::lowest_row_in_column(unsigned col) { column_strip & colstrip = m_columns[col]; - lp_assert(colstrip.size() > 0); + lp_assert(colstrip.live_size() > 0); unsigned ret = 0; - for (auto & t : colstrip) { - if (t.m_i > ret) { - ret = t.m_i; + for (const auto & t : colstrip) { + if (t.dead()) continue; + if (t.var() > ret) { + ret = t.var(); } } return ret; @@ -136,10 +146,10 @@ template void static_matrix::forget_last_colum template void static_matrix::remove_last_column(unsigned j) { column_strip & col = m_columns.back(); for (auto & it : col) { - auto & row = m_rows[it.m_i]; + auto & row = m_rows[it.var()]; unsigned offset = row.size() - 1; for (auto row_it = row.rbegin(); row_it != row.rend(); row_it ++) { - if (row_it->m_j == j) { + if (row_it->var() == j) { row.erase(row.begin() + offset); break; } @@ -150,24 +160,12 @@ template void static_matrix::remove_last_column(u m_vector_of_row_offsets.pop_back(); } - - - -template void static_matrix::set(unsigned row, unsigned col, T const & val) { - if (numeric_traits::is_zero(val)) return; - lp_assert(row < row_count() && col < column_count()); - auto & r = m_rows[row]; - unsigned offs_in_cols = static_cast(m_columns[col].size()); - m_columns[col].push_back(make_column_cell(row, static_cast(r.size()))); - r.push_back(make_row_cell(col, offs_in_cols, val)); -} - template std::set> static_matrix::get_domain() { std::set> ret; for (unsigned i = 0; i < m_rows.size(); i++) { for (auto &it : m_rows[i]) { - ret.insert(std::make_pair(i, it.m_j)); + ret.insert(std::make_pair(i, it.var())); } } return ret; @@ -179,7 +177,7 @@ template void static_matrix::copy_column_to_in for (auto & it : m_columns[j]) { const T& val = get_val(it); if (!is_zero(val)) - v.set_value(val, it.m_i); + v.set_value(val, it.var()); } } @@ -243,7 +241,7 @@ template void static_matrix::check_consistency std::unordered_map, T> by_rows; for (int i = 0; i < m_rows.size(); i++){ for (auto & t : m_rows[i]) { - std::pair p(i, t.m_j); + std::pair p(i, t.var()); lp_assert(by_rows.find(p) == by_rows.end()); by_rows[p] = t.get_val(); } @@ -251,7 +249,7 @@ template void static_matrix::check_consistency std::unordered_map, T> by_cols; for (int i = 0; i < m_columns.size(); i++){ for (auto & t : m_columns[i]) { - std::pair p(t.m_i, i); + std::pair p(t.var(), i); lp_assert(by_cols.find(p) == by_cols.end()); by_cols[p] = get_val(t); } @@ -266,51 +264,10 @@ template void static_matrix::check_consistency } #endif - -template void static_matrix::cross_out_row(unsigned k) { -#ifdef Z3DEBUG - check_consistency(); -#endif - cross_out_row_from_columns(k, m_rows[k]); - fix_row_indices_in_each_column_for_crossed_row(k); - m_rows.erase(m_rows.begin() + k); -#ifdef Z3DEBUG - regen_domain(); - check_consistency(); -#endif -} - - -template void static_matrix::fix_row_indices_in_each_column_for_crossed_row(unsigned k) { - for (unsigned j = 0; j < m_columns.size(); j++) { - auto & col = m_columns[j]; - for (int i = 0; i < col.size(); i++) { - if (col[i].m_i > k) { - col[i].m_i--; - } - } - } -} - -template void static_matrix::cross_out_row_from_columns(unsigned k, row_strip & row) { - for (auto & t : row) { - cross_out_row_from_column(t.m_j, k); - } -} - -template void static_matrix::cross_out_row_from_column(unsigned col, unsigned k) { - auto & s = m_columns[col]; - for (unsigned i = 0; i < s.size(); i++) { - if (s[i].m_i == k) { - s.erase(s.begin() + i); - break; - } - } -} - template T static_matrix::get_elem(unsigned i, unsigned j) const { // should not be used in efficient code !!!! for (auto & t : m_rows[i]) { - if (t.m_j == j) { + if (t.dead()) continue; + if (t.var() == j) { return t.get_val(); } } @@ -342,16 +299,22 @@ template bool static_matrix::is_correct() const { auto &r = m_rows[i]; std::unordered_set s; for (auto & rc : r) { - if (s.find(rc.m_j) != s.end()) { + if (rc.dead()) continue; + if (s.find(rc.var()) != s.end()) { return false; } - s.insert(rc.m_j); - if (rc.m_j >= m_columns.size()) + s.insert(rc.var()); + if (rc.var() >= m_columns.size()) return false; - if (rc.m_offset >= m_columns[rc.m_j].size()) + const auto& col = m_columns[rc.var()]; + if (col.cells_size() <= rc.offset()) return false; - if (rc.get_val() != get_val(m_columns[rc.m_j][rc.m_offset])) + const auto &cc = col[rc.offset()]; + if (cc.dead()) return false; + if (& m_rows[cc.var()][cc.offset()] != & rc) { + return false; + } if (is_zero(rc.get_val())) { return false; } @@ -363,51 +326,56 @@ template bool static_matrix::is_correct() const { auto & c = m_columns[j]; std::unordered_set s; for (auto & cc : c) { - if (s.find(cc.m_i) != s.end()) { + if (cc.dead()) + continue; + if (s.find(cc.var()) != s.end()) { return false; } - s.insert(cc.m_i); - if (cc.m_i >= m_rows.size()) + s.insert(cc.var()); + if (cc.var() >= m_rows.size()) return false; - if (cc.m_offset >= m_rows[cc.m_i].size()) + auto & rc = m_rows[cc.var()][cc.offset()]; + if (rc.dead()) return false; - if (get_val(cc) != m_rows[cc.m_i][cc.m_offset].get_val()) + if (&cc != &m_columns[rc.var()][rc.offset()]) return false; } } - + for (auto & row: m_rows) { + if (! row.is_correct()) + return false; + } + for (auto & col: m_columns) { + if (! col.is_correct()) + return false; + } return true; } template -void static_matrix::remove_element(vector> & row_vals, row_cell & row_el_iv) { - unsigned column_offset = row_el_iv.m_offset; - auto & column_vals = m_columns[row_el_iv.m_j]; - column_cell& cs = m_columns[row_el_iv.m_j][column_offset]; - unsigned row_offset = cs.m_offset; - if (column_offset != column_vals.size() - 1) { - auto & cc = column_vals[column_offset] = column_vals.back(); // copy from the tail - m_rows[cc.m_i][cc.m_offset].m_offset = column_offset; - } - - if (row_offset != row_vals.size() - 1) { - auto & rc = row_vals[row_offset] = row_vals.back(); // copy from the tail - m_columns[rc.m_j][rc.m_offset].m_offset = row_offset; - } - - column_vals.pop_back(); - row_vals.pop_back(); +void static_matrix::remove_element(row_cell & rc) { + lp_assert(rc.alive()); + unsigned j = rc.var(); + unsigned j_offset = rc.offset(); + auto & col = m_columns[j]; + column_cell & c = col[j_offset]; + unsigned i = c.var(); + unsigned i_offset = c.offset(); + col.delete_at(j_offset); + m_rows[i].delete_at(i_offset); } + template -void static_matrix::add_new_element(unsigned row, unsigned col, const T& val) { - auto & row_vals = m_rows[row]; - auto & col_vals = m_columns[col]; - unsigned row_el_offs = static_cast(row_vals.size()); - unsigned col_el_offs = static_cast(col_vals.size()); - row_vals.push_back(row_cell(col, col_el_offs, val)); - col_vals.push_back(column_cell(row, row_el_offs)); +void static_matrix::add_new_element(unsigned i, unsigned j, const T& val) { + auto & row = m_rows[i]; + auto & col = m_columns[j]; + unsigned offset_in_row, offset_in_col; + row_cell& rc = row.add_cell(j, val, offset_in_row); + column_cell& cc = col.add_cell(i, offset_in_col); + rc.offset() = offset_in_col; + cc.offset() = offset_in_row; } } From e9595eb2836a8878000d6b77c08c5b88c02a948a Mon Sep 17 00:00:00 2001 From: Lev Date: Sun, 29 Jul 2018 21:15:42 -0700 Subject: [PATCH 2/7] merge with z3prover Signed-off-by: Lev --- src/tactic/smtlogics/qflia_tactic.cpp | 4 --- src/util/lp/bound_propagator.cpp | 2 -- src/util/lp/lp_core_solver_base_def.h | 2 -- src/util/lp/lp_primal_core_solver.h | 24 +++++++-------- src/util/lp/static_matrix.h | 44 ++------------------------- 5 files changed, 13 insertions(+), 63 deletions(-) diff --git a/src/tactic/smtlogics/qflia_tactic.cpp b/src/tactic/smtlogics/qflia_tactic.cpp index ad5b660c3..eed4e4425 100644 --- a/src/tactic/smtlogics/qflia_tactic.cpp +++ b/src/tactic/smtlogics/qflia_tactic.cpp @@ -212,9 +212,6 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) { tactic * st = using_params(and_then(preamble_st, -#if 1 - mk_smt_tactic()), -#else or_else(mk_ilp_model_finder_tactic(m), mk_pb_tactic(m), and_then(fail_if_not(mk_is_quasi_pb_probe()), @@ -222,7 +219,6 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) { mk_fail_if_undecided_tactic()), mk_bounded_tactic(m), mk_smt_tactic())), -#endif main_p); diff --git a/src/util/lp/bound_propagator.cpp b/src/util/lp/bound_propagator.cpp index 199fef9b0..c4fa2aefa 100644 --- a/src/util/lp/bound_propagator.cpp +++ b/src/util/lp/bound_propagator.cpp @@ -16,8 +16,6 @@ const impq & bound_propagator::get_upper_bound(unsigned j) const { return m_lar_solver.m_mpq_lar_core_solver.m_r_upper_bounds()[j]; } void bound_propagator::try_add_bound(mpq v, unsigned j, bool is_low, bool coeff_before_j_is_pos, unsigned row_or_term_index, bool strict) { - TRACE("try_add_bound", - tout << "v = " << v << ", j = " << j << std::endl;); j = m_lar_solver.adjust_column_index_to_term_index(j); if (m_lar_solver.is_term(j)) { // lp treats terms as not having a free coefficient, restoring it below for the outside consumption diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index ddae7bf20..f9c1ae632 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -108,8 +108,6 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) { if (r.var() != j) m_d[r.var()] -= a * r.get_val(); } -ls - a = zero_of_type(); // zero the pivot column's m_d finally } diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index 2700fc168..d93174bc4 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -280,17 +280,14 @@ public: 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; + vector*> choices; + 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].m_cells.size(); k++) { const row_cell& rc = this->m_A.m_rows[i].m_cells[k]; if (rc.dead()) continue; - const row_cell& rc = this->m_A.m_rows[i].m_cells[k]; - if (rc.dead()) continue; ->>>>>>> e6c612f... trying the new scheme in static_matrix : in progress unsigned j = rc.var(); if (j == bj) continue; @@ -305,23 +302,24 @@ public: if (damage < num_of_non_free_basics) { num_of_non_free_basics = damage; len = this->m_A.m_columns[j].live_size(); - choice = k; - nchoices = 1; + choices.clear(); + choices.push_back(&rc); } else if (damage == num_of_non_free_basics && - this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % (++nchoices))) { - choice = k; + this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % 2)) { + choices.push_back(&rc); len = this->m_A.m_columns[j].live_size(); } } - if (choice == -1) { + if (choices.size() == 0) { m_inf_row_index_for_tableau = i; return -1; } - const row_cell& rc = this->m_A.m_rows[i].m_cells[choice]; - a_ent = rc.coeff(); - return rc.var(); + const row_cell* rc = choices.size() == 1? choices[0] : + choices[this->m_settings.random_next() % choices.size()]; + a_ent = rc->coeff(); + return rc->var(); } static X positive_infinity() { return convert_struct::convert(std::numeric_limits::max()); diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index 036f628b7..c75fa4ce8 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -64,9 +64,6 @@ public: bool alive() const { return !dead(); } bool dead() const { return m_i == static_cast(-1); } void set_dead() { m_i = -1;} - - - }; template @@ -79,6 +76,7 @@ public: } row_cell(unsigned j, T const & val) : m_j(j), m_value(val) { } + const T & get_val() const { lp_assert(alive()); return m_value; @@ -126,8 +124,6 @@ public: bool alive() const { return !dead(); } bool dead() const { return m_j == static_cast(-1); } void set_dead() { m_j = static_cast(-1); } - - }; template @@ -252,10 +248,6 @@ public: return m_live_size; } - unsigned cells_size() const { - return m_cells.size(); - } - unsigned cells_size() const { return m_cells.size(); } @@ -304,7 +296,7 @@ public: } m_cells.pop_back(); } - + bool is_correct() const { std::set d0; std::set d1; unsigned alive = 0; @@ -341,38 +333,6 @@ public: lp_assert(is_correct()); } - void swap_with_head_cell(unsigned i) { - lp_assert(i > 0); - lp_assert(m_cells[i].alive()); - column_cell head_copy = m_cells[0]; - if (head_copy.dead()) { - if (m_first_dead == 0) { - m_first_dead = i; - } else { - column_cell * c = &m_cells[m_first_dead]; - for (; c->next_dead_index() != 0; c = &m_cells[c->next_dead_index()]); - lp_assert(c->next_dead_index() == 0); - c->next_dead_index() = i; - } - } - m_cells[0] = m_cells[i]; - m_cells[i] = head_copy; - lp_assert(is_correct()); - } - - column_cell & add_cell(unsigned i, unsigned & index) { - if (m_first_dead != -1) { - auto & ret = m_cells[index = m_first_dead]; - m_first_dead = ret.next_dead_index(); - m_live_size++; - ret.var() = i; - return ret; - } - lp_assert(m_live_size == m_cells.size()); - index = m_live_size++; - m_cells.push_back(column_cell(i)); - return m_cells.back(); - } column_cell & add_cell(unsigned i, unsigned & index) { if (m_first_dead != -1) { auto & ret = m_cells[index = m_first_dead]; From 2de27ae3afb727de53b135ec542d51b9205c56a3 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Sun, 29 Jul 2018 22:33:19 -0700 Subject: [PATCH 3/7] uniform choice of a beneficial column Signed-off-by: Lev Nachmanson --- src/util/lp/lp_primal_core_solver.h | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index d93174bc4..faf9074f6 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -280,7 +280,8 @@ public: 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 - vector*> choices; + int choice = -1; + int nchoices = 0; unsigned num_of_non_free_basics = 1000000; unsigned len = 100000000; unsigned bj = this->m_basis[i]; @@ -302,24 +303,23 @@ public: if (damage < num_of_non_free_basics) { num_of_non_free_basics = damage; len = this->m_A.m_columns[j].live_size(); - choices.clear(); - choices.push_back(&rc); + choice = k; + nchoices = 1; } else if (damage == num_of_non_free_basics && - this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % 2)) { - choices.push_back(&rc); + this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % (++nchoices))) { + choice = k; len = this->m_A.m_columns[j].live_size(); } } - if (choices.size() == 0) { + if (choice == -1) { m_inf_row_index_for_tableau = i; return -1; } - const row_cell* rc = choices.size() == 1? choices[0] : - choices[this->m_settings.random_next() % choices.size()]; - a_ent = rc->coeff(); - return rc->var(); + const row_cell& rc = this->m_A.m_rows[i].m_cells[choice]; + a_ent = rc.coeff(); + return rc.var(); } static X positive_infinity() { return convert_struct::convert(std::numeric_limits::max()); From 9cb713879ef367d44e65ed6dc7089f172ded4348 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 30 Jul 2018 09:56:39 -0700 Subject: [PATCH 4/7] fix the build Signed-off-by: Lev Nachmanson --- src/util/lp/lp_core_solver_base.h | 8 -------- 1 file changed, 8 deletions(-) diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 35a989c07..5bde78ca8 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -648,14 +648,6 @@ public: return true; } - int find_pivot_index_in_row(unsigned i, const vector & col) const { - for (const auto & c: col) { - if (c.m_i == i) - return c.m_offset; - } - return -1; - } - void transpose_rows_tableau(unsigned i, unsigned ii); void pivot_to_reduced_costs_tableau(unsigned i, unsigned j); From 181bb60e363e5aaebd87d9286465da8da2ee53a1 Mon Sep 17 00:00:00 2001 From: Lev Date: Mon, 30 Jul 2018 12:54:53 -0700 Subject: [PATCH 5/7] remove some lp_asserts Signed-off-by: Lev --- src/util/lp/lp_core_solver_base.h | 1 - src/util/lp/lp_core_solver_base_def.h | 1 - src/util/lp/static_matrix.h | 2 -- src/util/lp/static_matrix_def.h | 4 ++-- 4 files changed, 2 insertions(+), 6 deletions(-) diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 5bde78ca8..fbb08edfe 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -209,7 +209,6 @@ public: bool need_to_pivot_to_basis_tableau() const { - lp_assert(m_A.is_correct()); unsigned m = m_A.row_count(); for (unsigned i = 0; i < m; i++) { unsigned bj = m_basis[i]; diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index f9c1ae632..8cfc39170 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -621,7 +621,6 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) { } template bool lp_core_solver_base:: pivot_column_tableau(unsigned j, unsigned piv_row_index) { - lp_assert(m_A.is_correct()); m_A.compress_row_if_needed(piv_row_index); if (!divide_row_by_pivot(piv_row_index, j)) return false; diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index c75fa4ce8..e3ea3c1b0 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -790,12 +790,10 @@ public: void compress_row_if_needed(unsigned i) { compress_cells(m_rows[i], m_columns); - lp_assert(is_correct()); } void compress_column_if_needed(unsigned j) { compress_cells(m_columns[j], m_rows); - lp_assert(is_correct()); } ref_row operator[](unsigned i) const { return ref_row(*this, i);} diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index 698363289..ae0e22540 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -44,7 +44,6 @@ template void static_matrix::scan_row_ii_to_offse template bool static_matrix::pivot_row_to_row_given_cell(unsigned i, column_cell & c, unsigned pivot_col) { - lp_assert(is_correct()); unsigned ii = c.var(); lp_assert(i < row_count() && ii < column_count() && i != ii); T alpha = -get_val(c); @@ -83,7 +82,6 @@ template bool static_matrix::pivot_row_to_row_giv if (is_zero(c.coeff())) remove_element(c); } - lp_assert(is_correct()); return !rowii.empty(); } @@ -295,6 +293,8 @@ template T static_matrix::get_row_balance(unsi } template bool static_matrix::is_correct() const { + if (m_rows.size() > 100 || m_columns.size() > 100) + return true; for (unsigned i = 0; i < m_rows.size(); i++) { auto &r = m_rows[i]; std::unordered_set s; From 0ee68220e1ff210e663434cc3559b9486135188a Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Mon, 30 Jul 2018 14:34:03 -0700 Subject: [PATCH 6/7] use CASSERT instead of lp_assert for static_matrix Signed-off-by: Lev Nachmanson --- src/util/lp/lar_solver.cpp | 2 +- src/util/lp/lp_core_solver_base.h | 2 +- src/util/lp/lp_core_solver_base_def.h | 6 +++--- src/util/lp/static_matrix.h | 8 ++++---- src/util/lp/static_matrix_def.h | 2 -- 5 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 5408bd922..0d8e7f70a 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -1373,7 +1373,7 @@ void lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { lp_assert(A_r().m_columns[j].live_size() == 0); A_r().m_rows.pop_back(); A_r().m_columns.pop_back(); - lp_assert(A_r().is_correct()); + CASSERT("check_static_matrix", A_r().is_correct()); slv.m_b.pop_back(); } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index fbb08edfe..00ba437d1 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -230,7 +230,7 @@ public: bool reduced_costs_are_correct_tableau() const { if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) return true; - lp_assert(m_A.is_correct()); + CASSERT("check_static_matrix", m_A.is_correct()); if (m_using_infeas_costs) { if (infeasibility_costs_are_correct() == false) { return false; diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index 8cfc39170..6a9af8a26 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -616,7 +616,7 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) { } } coeff = one_of_type(); - lp_assert(m_A.is_correct()); + CASSERT("check_static_matrix", m_A.is_correct()); return true; } template bool lp_core_solver_base:: @@ -639,7 +639,7 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { if (pivot_col_cell_index != 0) m_A.swap_with_head_cell(j, pivot_col_cell_index); - lp_assert(m_A.is_correct()); + CASSERT("check_static_matrix", m_A.is_correct()); while (column.live_size() > 1) { auto & c = column.back(); if (c.dead()) { @@ -655,7 +655,7 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { } m_A.compress_column_if_needed(j); lp_assert(column.live_size() == 1); - lp_assert(m_A.is_correct()); + CASSERT("check_static_matrix", m_A.is_correct()); if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs) pivot_to_reduced_costs_tableau(piv_row_index, j); diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index e3ea3c1b0..cb6af1379 100644 --- a/src/util/lp/static_matrix.h +++ b/src/util/lp/static_matrix.h @@ -184,7 +184,7 @@ public: c.next_dead_index() = m_first_dead; m_first_dead = j; } - lp_assert(is_correct()); + CASSERT("check_static_matrix", is_correct()); } bool is_correct() const { @@ -264,7 +264,7 @@ public: c.next_dead_index() = m_first_dead; m_first_dead = j; } - lp_assert(is_correct()); + CASSERT("check_static_matrix", is_correct()); } const column_cell& operator[] (unsigned i) const { return m_cells[i];} @@ -330,7 +330,7 @@ public: } m_cells[0] = m_cells[i]; m_cells[i] = head_copy; - lp_assert(is_correct()); + CASSERT("check_static_matrix", is_correct()); } column_cell & add_cell(unsigned i, unsigned & index) { @@ -580,7 +580,7 @@ public: m_columns.pop_back(); // delete the last column m_stack.pop(); } - lp_assert(is_correct()); + CASSERT("check_static_matrix", is_correct()); } void multiply_row(unsigned row, T const & alpha) { diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index ae0e22540..43c91cefa 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -293,8 +293,6 @@ template T static_matrix::get_row_balance(unsi } template bool static_matrix::is_correct() const { - if (m_rows.size() > 100 || m_columns.size() > 100) - return true; for (unsigned i = 0; i < m_rows.size(); i++) { auto &r = m_rows[i]; std::unordered_set s; From 3d274c2e6f58ec4500c0c6f355e007f28d0a03a6 Mon Sep 17 00:00:00 2001 From: Lev Date: Mon, 30 Jul 2018 15:55:06 -0700 Subject: [PATCH 7/7] use CASSERT for hnf Signed-off-by: Lev --- src/util/lp/hnf.h | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/util/lp/hnf.h b/src/util/lp/hnf.h index 3cdeac466..8edce39e0 100644 --- a/src/util/lp/hnf.h +++ b/src/util/lp/hnf.h @@ -317,7 +317,6 @@ class hnf { } void handle_column_ij_in_row_i(unsigned i, unsigned j) { - lp_assert(is_correct_modulo()); const mpq& aii = m_H[i][i]; const mpq& aij = m_H[i][j]; mpq p,q,r; @@ -342,7 +341,7 @@ class hnf { // from the left multiply_U_reverse_from_left_by(i, j, aii_over_r, aij_over_r, -q, p); - lp_assert(is_correct_modulo()); + CASSERT("hnf_test", is_correct_modulo()); } @@ -402,8 +401,6 @@ class hnf { } void process_row(unsigned i) { - - lp_assert(is_correct_modulo()); for (unsigned j = i + 1; j < m_n; j++) { process_row_column(i, j); } @@ -414,7 +411,7 @@ class hnf { if (is_neg(m_H[i][i])) switch_sign_for_column(i); work_on_columns_less_than_i_in_the_triangle(i); - lp_assert(is_correct_modulo()); + CASSERT("hnf_test", is_correct_modulo()); } void calculate() { @@ -604,7 +601,7 @@ public: #ifdef Z3DEBUG prepare_U_and_U_reverse(); calculate(); - lp_assert(is_correct_final()); + CASSERT("hnf_test", is_correct_final()); #endif calculate_by_modulo(); #ifdef Z3DEBUG