From 0a5141780431b8f4af8e4c44ca4da5aa864645a1 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Tue, 31 Jul 2018 22:51:27 -0700 Subject: [PATCH 1/4] unroll static_matrix to avoid dead cells Signed-off-by: Lev Nachmanson --- src/util/lp/bound_analyzer_on_row.h | 8 - src/util/lp/int_solver.cpp | 11 - src/util/lp/lar_core_solver.h | 5 +- src/util/lp/lar_core_solver_def.h | 3 +- src/util/lp/lar_solver.cpp | 29 +- src/util/lp/lp_core_solver_base.h | 7 +- src/util/lp/lp_core_solver_base_def.h | 44 +- src/util/lp/lp_primal_core_solver.h | 22 +- src/util/lp/lp_primal_core_solver_def.h | 2 +- .../lp/lp_primal_core_solver_tableau_def.h | 11 +- src/util/lp/random_updater_def.h | 3 +- src/util/lp/static_matrix.cpp | 9 +- src/util/lp/static_matrix.h | 502 +++--------------- src/util/lp/static_matrix_def.h | 168 +++--- 14 files changed, 231 insertions(+), 593 deletions(-) diff --git a/src/util/lp/bound_analyzer_on_row.h b/src/util/lp/bound_analyzer_on_row.h index 3d047cd86..549c8e5ce 100644 --- a/src/util/lp/bound_analyzer_on_row.h +++ b/src/util/lp/bound_analyzer_on_row.h @@ -59,8 +59,6 @@ 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()); @@ -170,7 +168,6 @@ 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) @@ -179,7 +176,6 @@ 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); @@ -197,7 +193,6 @@ 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) @@ -205,7 +200,6 @@ 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); @@ -228,7 +222,6 @@ 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(); @@ -258,7 +251,6 @@ 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/int_solver.cpp b/src/util/lp/int_solver.cpp index 78206dc53..a77c202a0 100644 --- a/src/util/lp/int_solver.cpp +++ b/src/util/lp/int_solver.cpp @@ -100,7 +100,6 @@ 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)) @@ -312,7 +311,6 @@ 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; @@ -327,7 +325,6 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip & row mpq f_0 = int_solver::fractional_part(get_value(inf_col)); mpq one_min_f_0 = 1 - f_0; for (const auto & p : row) { - if (p.dead()) continue; x_j = p.var(); if (x_j == inf_col) continue; @@ -356,7 +353,6 @@ 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); @@ -793,7 +789,6 @@ 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; @@ -807,7 +802,6 @@ 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)) { @@ -873,7 +867,6 @@ 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()); @@ -899,7 +892,6 @@ 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)) @@ -1031,7 +1023,6 @@ 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(); @@ -1162,14 +1153,12 @@ 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 (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 (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 4223f5964..904550339 100644 --- a/src/util/lp/lar_core_solver.h +++ b/src/util/lp/lar_core_solver.h @@ -583,7 +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].live_size() == 1); + lp_assert(m_r_solver.m_A.m_columns[j].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; @@ -632,8 +632,7 @@ 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.m_cells) { - if (c.dead()) continue; + for (row_cell & c : row) { 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 c626362fc..b945b0e51 100644 --- a/src/util/lp/lar_core_solver_def.h +++ b/src/util/lp/lar_core_solver_def.h @@ -226,8 +226,7 @@ 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_cells) { - if (rc.dead()) continue; + 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.var())); } } diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 0d8e7f70a..78bd45e94 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -162,7 +162,7 @@ void lar_solver::analyze_new_bounds_on_row_tableau( unsigned row_index, bound_propagator & bp ) { - if (A_r().m_rows[row_index].live_size() > settings().max_row_length_for_bound_propagation) + if (A_r().m_rows[row_index].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], @@ -215,7 +215,6 @@ 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]) { - if (r.dead()) continue; unsigned j = r.var(); if (j == bound_j) continue; mpq const& a = r.get_val(); @@ -429,7 +428,6 @@ void lar_solver::set_costs_to_zero(const lar_term& term) { jset.insert(j); else { for (const auto & rc : A_r().m_rows[i]) - if (rc.alive()) jset.insert(rc.var()); } } @@ -640,8 +638,7 @@ 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]) - if (rc.alive()) - m_rows_with_changed_bounds.insert(rc.var()); + m_rows_with_changed_bounds.insert(rc.var()); } bool lar_solver::use_tableau() const { return m_settings.use_tableau(); } @@ -670,7 +667,6 @@ 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]) { - if (c.dead()) continue; r += c.coeff() * m_mpq_lar_core_solver.m_r_x[c.var()]; } return is_zero(r); @@ -695,7 +691,6 @@ 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]) { - 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); @@ -871,7 +866,7 @@ 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].live_size() == 0); + lp_assert(A.m_rows[last_row].size() == 0); for (auto & t : ls->m_coeffs) { lp_assert(!is_zero(t.second)); var_index j = t.first; @@ -1332,9 +1327,8 @@ 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.cells_size(); k-- > 0;){ + for (unsigned k = last_column.size(); k-- > 0;){ auto & cc = last_column[k]; - if (cc.dead()) continue; if (cc.var() == i) return; non_zero_column_cell_index = k; @@ -1357,20 +1351,15 @@ 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.cells_size(); k-- > 0;) { + for (unsigned k = last_row.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.var()] += cost_j*rc.get_val(); } - - A_r().remove_element(rc); + A_r().remove_element(last_row, rc); } - lp_assert(last_row.live_size() == 0); - lp_assert(A_r().m_columns[j].live_size() == 0); + lp_assert(last_row.size() == 0); + lp_assert(A_r().m_columns[j].size() == 0); A_r().m_rows.pop_back(); A_r().m_columns.pop_back(); CASSERT("check_static_matrix", A_r().is_correct()); @@ -1379,7 +1368,7 @@ void lar_solver::remove_last_row_and_column_from_tableau(unsigned j) { void lar_solver::remove_last_column_from_A() { // the last column has to be empty - lp_assert(A_r().m_columns.back().live_size() == 0); + lp_assert(A_r().m_columns.back().size() == 0); A_r().m_columns.pop_back(); } diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h index 00ba437d1..41b6fe31d 100644 --- a/src/util/lp/lp_core_solver_base.h +++ b/src/util/lp/lp_core_solver_base.h @@ -212,12 +212,10 @@ 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].live_size() > 0); - if (m_A.m_columns[bj].live_size() > 1) + lp_assert(m_A.m_columns[bj].size() > 0); + if (m_A.m_columns[bj].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 @@ -246,7 +244,6 @@ public: } else { auto d = m_costs[j]; 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) { diff --git a/src/util/lp/lp_core_solver_base_def.h b/src/util/lp/lp_core_solver_base_def.h index 6a9af8a26..c3a0a0a00 100644 --- a/src/util/lp/lp_core_solver_base_def.h +++ b/src/util/lp/lp_core_solver_base_def.h @@ -104,7 +104,6 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) { if (is_zero(a)) return; 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(); } @@ -309,8 +308,7 @@ 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].m_cells) { - if (c.dead()) continue; + for (auto & c : m_A.m_rows[i]) { 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); @@ -591,10 +589,9 @@ 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.cells_size(); + unsigned size = row.size(); for (unsigned j = 0; j < size; j++) { auto & c = row[j]; - if (c.dead()) continue; if (c.var() == pivot_col) { pivot_index = static_cast(j); break; @@ -610,7 +607,6 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) { this->m_b[pivot_row] /= coeff; for (unsigned j = 0; j < size; j++) { auto & c = row[j]; - if (c.dead()) continue; if (c.var() != pivot_col) { c.coeff() /= coeff; } @@ -621,13 +617,11 @@ 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) { - m_A.compress_row_if_needed(piv_row_index); - if (!divide_row_by_pivot(piv_row_index, j)) + 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.cells_size(); k++) { - if (column[k].dead()) continue; + for (unsigned k = 0; k < column.size(); k++) { if (column[k].var() == piv_row_index) { pivot_col_cell_index = k; break; @@ -636,16 +630,18 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { if (pivot_col_cell_index < 0) return false; - if (pivot_col_cell_index != 0) - m_A.swap_with_head_cell(j, pivot_col_cell_index); + 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; - CASSERT("check_static_matrix", m_A.is_correct()); - while (column.live_size() > 1) { + m_A.m_rows[piv_row_index][column[0].m_offset].m_offset = 0; + m_A.m_rows[c.var()][c.m_offset].m_offset = pivot_col_cell_index; + } + while (column.size() > 1) { auto & c = column.back(); - 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; @@ -653,9 +649,6 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) { if (m_pivoted_rows!= nullptr) m_pivoted_rows->insert(c.var()); } - m_A.compress_column_if_needed(j); - lp_assert(column.live_size() == 1); - 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); @@ -769,8 +762,7 @@ fill_reduced_costs_from_m_y_by_rows() { while (i--) { const T & y = m_y[i]; if (is_zero(y)) continue; - for (row_cell & c : m_A.m_rows[i].m_cells) { - if (c.dead()) continue; + for (row_cell & c : m_A.m_rows[i]) { j = c.var(); if (m_basis_heading[j] < 0) { m_d[j] -= y * c.get_val(); @@ -969,8 +961,7 @@ 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].m_cells) { - if (c.dead()) continue; + for (auto &c : m_A.m_rows[i]) { j = c.var(); if (j == basic_j) continue; @@ -985,8 +976,7 @@ 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].m_cells) { - if (c.dead()) continue; + for (auto &c : m_A.m_rows[i]) { if (c.var() == basic_j) continue; if (pivot_column_general(c.var(), basic_j, w)) diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h index faf9074f6..db92edb05 100644 --- a/src/util/lp/lp_primal_core_solver.h +++ b/src/util/lp/lp_primal_core_solver.h @@ -240,7 +240,6 @@ 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 (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++; @@ -254,7 +253,6 @@ 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.dead()) continue; if (rc.var() == bj) continue; if (bj_needs_to_grow) { @@ -286,9 +284,8 @@ public: 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; + 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; @@ -302,13 +299,13 @@ 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].live_size(); + 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].live_size() <= len && (this->m_settings.random_next() % (++nchoices))) { + this->m_A.m_columns[j].size() <= len && (this->m_settings.random_next() % (++nchoices))) { choice = k; - len = this->m_A.m_columns[j].live_size(); + len = this->m_A.m_columns[j].size(); } } @@ -317,7 +314,7 @@ public: m_inf_row_index_for_tableau = i; return -1; } - const row_cell& rc = this->m_A.m_rows[i].m_cells[choice]; + const row_cell& rc = this->m_A.m_rows[i][choice]; a_ent = rc.coeff(); return rc.var(); } @@ -830,11 +827,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].live_size(); + this->m_columns_nz[i] = this->m_A.m_columns[i].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].live_size(); + this->m_rows_nz[i] = this->m_A.m_rows[i].size(); } } @@ -864,7 +861,7 @@ public: unsigned solve_with_tableau(); bool basis_column_is_set_correctly(unsigned j) const { - return this->m_A.m_columns[j].live_size() == 1; + return this->m_A.m_columns[j].size() == 1; } @@ -884,7 +881,6 @@ 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]) { - if (rc.dead()) continue; unsigned k = rc.var(); if (k == j) continue; diff --git a/src/util/lp/lp_primal_core_solver_def.h b/src/util/lp/lp_primal_core_solver_def.h index eb4a9ce1d..1e9edbd31 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].live_size() + 1)) + this->m_column_norms[j] = T(static_cast(this->m_A.m_columns[j].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 be28bf36f..8f9cf7ff1 100644 --- a/src/util/lp/lp_primal_core_solver_tableau_def.h +++ b/src/util/lp/lp_primal_core_solver_tableau_def.h @@ -266,10 +266,9 @@ 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.cells_size(); + unsigned col_size = col.size(); for (;k < col_size && unlimited; k++) { const column_cell & c = col[k]; - if (c.dead()) continue; unsigned i = c.var(); const T & ed = this->m_A.get_val(c); lp_assert(!numeric_traits::is_zero(ed)); @@ -277,7 +276,7 @@ template int lp_primal_core_solver::find_leaving_ 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].live_size(); + row_min_nz = this->m_A.m_rows[i].size(); } } if (unlimited) { @@ -289,7 +288,6 @@ template int lp_primal_core_solver::find_leaving_ X ratio; for (;k < col_size; k++) { const column_cell & c = col[k]; - if (c.dead()) continue; unsigned i = c.var(); const T & ed = this->m_A.get_val(c); lp_assert(!numeric_traits::is_zero(ed)); @@ -297,7 +295,7 @@ template int lp_primal_core_solver::find_leaving_ 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].live_size(); + unsigned i_nz = this->m_A.m_rows[i].size(); if (ratio < t) { t = ratio; m_leaving_candidates.clear(); @@ -306,7 +304,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].live_size(); + row_min_nz = this->m_A.m_rows[i].size(); } else if (ratio == t && i_nz == row_min_nz) { m_leaving_candidates.push_back(j); } @@ -361,7 +359,6 @@ 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]) { - 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)); } diff --git a/src/util/lp/random_updater_def.h b/src/util/lp/random_updater_def.h index 06cc01163..c972cc175 100644 --- a/src/util/lp/random_updater_def.h +++ b/src/util/lp/random_updater_def.h @@ -83,8 +83,7 @@ void random_updater::add_column_to_sets(unsigned j) { add_value(m_lar_solver.get_core_solver().m_r_x[j]); } else { unsigned row = m_lar_solver.get_core_solver().m_r_heading[j]; - for (auto & row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row].m_cells) { - if (row_c.dead()) continue; + for (auto & row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row]) { unsigned cj = row_c.var(); if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) { m_var_set.insert(cj); diff --git a/src/util/lp/static_matrix.cpp b/src/util/lp/static_matrix.cpp index d0448b6c9..9d12087bb 100644 --- a/src/util/lp/static_matrix.cpp +++ b/src/util/lp/static_matrix.cpp @@ -30,9 +30,6 @@ Revision History: #include "util/lp/lar_solver.h" namespace lp { template void static_matrix::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; @@ -50,6 +47,7 @@ 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); @@ -65,6 +63,7 @@ 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 @@ -73,11 +72,13 @@ 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(lp::row_cell&); +template void lp::static_matrix >::remove_element(vector, true, unsigned int>&, lp::row_cell&); } diff --git a/src/util/lp/static_matrix.h b/src/util/lp/static_matrix.h index cb6af1379..bb6cdefc8 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,347 +29,31 @@ #include namespace lp { -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 -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: +struct row_cell { + unsigned m_j; // points to the column + unsigned m_offset; // offset in column + T m_value; row_cell(unsigned j, unsigned offset, T const & val) : m_j(j), m_offset(offset), m_value(val) { } - row_cell(unsigned j, T const & val) : m_j(j), m_value(val) { + row_cell(unsigned j, unsigned offset) : m_j(j), m_offset(offset) { } - - 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); } + const T & get_val() const { return m_value;} + T & get_val() { return m_value;} + const T & coeff() const { return m_value;} + T & coeff() { return m_value;} + unsigned var() const { return m_j;} + unsigned & var() { return m_j;} + unsigned offset() const { return m_offset;} + unsigned & offset() { return m_offset;} + }; +struct empty_struct {}; +typedef row_cell column_cell; + template -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; - } - CASSERT("check_static_matrix", 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(); - } - - 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; - } - CASSERT("check_static_matrix", 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(); - } - 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; - } - - 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; - CASSERT("check_static_matrix", 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(); - } - 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(); -} - +using row_strip = vector>; // each assignment for this matrix should be issued only once!!! template @@ -383,13 +67,13 @@ class static_matrix unsigned m_n; dim(unsigned m, unsigned n) :m_m(m), m_n(n) {} }; - // fields - std::stack m_stack; + std::stack m_stack; public: - vector m_vector_of_row_offsets; - indexed_vector m_work_vector; - vector> m_rows; - vector m_columns; + typedef vector column_strip; + 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; @@ -397,9 +81,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.add_new_element( m_row, m_col, v); return *this; } + ref & operator=(T const & v) { m_matrix.set( m_row, m_col, v); 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; } + ref operator=(ref & v) { m_matrix.set(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); } }; @@ -415,7 +99,7 @@ public: public: const T & get_val(const column_cell & c) const { - return m_rows[c.var()][c.offset()].get_val(); + return m_rows[c.var()][c.m_offset].get_val(); } column_cell & get_column_cell(const row_cell &rc) { @@ -424,7 +108,7 @@ public: void init_row_columns(unsigned m, unsigned n); - // constructor with no parameters + // constructor with no parameters static_matrix() {} // constructor @@ -459,11 +143,11 @@ public: void remove_last_column(unsigned j); - void remove_element(row_cell & elem_to_remove); + void remove_element(vector> & row, row_cell & elem_to_remove); void multiply_column(unsigned column, T const & alpha) { for (auto & t : m_columns[column]) { - auto & r = m_rows[t.var()].m_cells[t.offset()]; + auto & r = m_rows[t.var()][t.offset()]; r.coeff() *= alpha; } } @@ -482,6 +166,8 @@ 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(); @@ -490,8 +176,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 & c : m_columns[j]) { - v[c.var()] += a * get_val(c); + for (const auto & it : m_columns[j]) { + v[it.var()] += a * get_val(it); } } @@ -504,6 +190,9 @@ public: void check_consistency(); #endif + + void cross_out_row(unsigned k); + // void fix_row_indices_in_each_column_for_crossed_row(unsigned k); @@ -514,9 +203,9 @@ public: T get_elem(unsigned i, unsigned j) const; - unsigned number_of_non_zeroes_in_column(unsigned j) const { return m_columns[j].live_size(); } + unsigned number_of_non_zeroes_in_column(unsigned j) const { return m_columns[j].size(); } - unsigned number_of_non_zeroes_in_row(unsigned i) const { return m_rows[i].live_size(); } + unsigned number_of_non_zeroes_in_row(unsigned i) const { return m_rows[i].size(); } unsigned number_of_non_zeroes() const { unsigned ret = 0; @@ -549,12 +238,12 @@ public: m_stack.push(d); } - void pop_row_columns(const row_strip & row) { + void pop_row_columns(const vector> & row) { for (auto & c : row) { - if (c.dead()) { - continue; - } - m_columns[c.var()].delete_at(c.offset()); + unsigned j = c.m_j; + auto & col = m_columns[j]; + lp_assert(col[col.size() - 1].var() == m_rows.size() -1 ); // todo : start here!!!! + col.pop_back(); } } @@ -580,13 +269,12 @@ public: m_columns.pop_back(); // delete the last column m_stack.pop(); } - CASSERT("check_static_matrix", is_correct()); + lp_assert(is_correct()); } void multiply_row(unsigned row, T const & alpha) { - for (auto & t : m_rows[row].m_cells) { - if (t.alive()) - t.coeff() *= alpha; + for (auto & t : m_rows[row]) { + t.m_value *= alpha; } } @@ -599,8 +287,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 & c : m_columns[j]) { - ret += y[c.var()] * get_val(c); + for (auto & it : m_columns[j]) { + ret += y[it.var()] * get_val(it); // get_value_of_column_cell(it); } return ret; } @@ -614,20 +302,18 @@ public: m_rows[i] = m_rows[ii]; m_rows[ii] = t; // now fix the columns - for (const auto & rc : m_rows[i]) { - if (rc.dead()) continue; - column_cell & cc = m_columns[rc.var()][rc.offset()]; + for (auto & rc : m_rows[i]) { + column_cell & cc = m_columns[rc.m_j][rc.m_offset]; lp_assert(cc.var() == ii); cc.var() = i; } - for (const auto & rc : m_rows[ii]) { - if (rc.dead()) continue; - column_cell & cc = m_columns[rc.var()][rc.offset()]; + for (auto & rc : m_rows[ii]) { + column_cell & cc = m_columns[rc.m_j][rc.m_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) @@ -637,18 +323,17 @@ public: return; for (const auto & c : m_rows[row_index]) { - if (c.dead()) continue; - if (c.var() == j) { + if (c.m_j == j) { continue; } - T & wv = m_work_vector.m_data[c.var()]; + T & wv = m_work_vector.m_data[c.m_j]; bool was_zero = is_zero(wv); - wv -= alpha * c.coeff(); + wv -= alpha * c.m_value; if (was_zero) - m_work_vector.m_index.push_back(c.var()); + m_work_vector.m_index.push_back(c.m_j); else { if (is_zero(wv)) { - m_work_vector.erase_from_index(c.var()); + m_work_vector.erase_from_index(c.m_j); } } } @@ -666,7 +351,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()); @@ -682,10 +367,10 @@ public: unsigned last_row = row_count() - 1; for (unsigned j : m_work_vector.m_index) { - add_new_element(last_row, j, m_work_vector.m_data[j]); + set (last_row, j, m_work_vector.m_data[j]); } lp_assert(column_count() > 0); - add_new_element(last_row, column_count() - 1, one_of_type()); + set(last_row, column_count() - 1, one_of_type()); } void copy_column_to_vector (unsigned j, vector & v) const { @@ -702,8 +387,7 @@ public: L ret = zero_of_type(); lp_assert(row < m_rows.size()); for (auto & it : m_rows[row]) { - if (it.dead()) continue; - ret += w[it.var()] * it.coeff(); + ret += w[it.m_j] * it.get_val(); } return ret; } @@ -715,8 +399,8 @@ public: column_cell_plus(const column_cell & c, const static_matrix& A) : m_c(c), m_A(A) {} 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(); } + const T & coeff() const { return m_A.m_rows[var()][m_c.m_offset].get_val(); } + }; struct column_container { @@ -728,7 +412,7 @@ public: // fields const column_cell *m_c; const static_matrix& m_A; - const column_cell *m_end; + //typedefs @@ -742,19 +426,13 @@ 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; } @@ -772,30 +450,8 @@ 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); - } - - void compress_column_if_needed(unsigned j) { - compress_cells(m_columns[j], m_rows); - } - 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 43c91cefa..1dcaa34f8 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -36,10 +36,8 @@ 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.cells_size(); j++) { - if (rvals[j].dead()) continue; - m_vector_of_row_offsets[rvals[j].var()] = j; - } + for (unsigned j = 0; j < rvals.size(); j++) + m_vector_of_row_offsets[rvals[j].var()] = j; } @@ -48,39 +46,33 @@ template bool static_matrix::pivot_row_to_row_giv 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[c.offset()]); + remove_element(rowii, rowii[c.m_offset]); scan_row_ii_to_offset_vector(rowii); - unsigned prev_size_ii = rowii.cells_size(); + unsigned prev_size_ii = rowii.size(); // run over the pivot row and update row ii - for (const auto & iv : m_rows[i].m_cells) { - if (iv.dead()) continue; + for (const auto & iv : m_rows[i]) { unsigned j = iv.var(); if (j == pivot_col) continue; - lp_assert(!is_zero(iv.coeff())); - T alv = alpha * iv.coeff(); + T alv = alpha * iv.m_value; + lp_assert(!is_zero(iv.m_value)); 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].coeff() += alv; + rowii[j_offs].m_value += alv; } } // clean the work vector for (unsigned k = 0; k < prev_size_ii; k++) { - auto & c = rowii[k]; - if (c.dead()) continue; - m_vector_of_row_offsets[c.var()] = -1; + m_vector_of_row_offsets[rowii[k].var()] = -1; } + // remove zeroes - for (unsigned k = rowii.cells_size(); k-- > 0; ) { - auto & c = rowii[k]; - if (c.dead()) - continue; - if (is_zero(c.coeff())) - remove_element(c); + for (unsigned k = rowii.size(); k-- > 0; ) { + if (is_zero(rowii[k].m_value)) + remove_element(rowii, rowii[k]); } return !rowii.empty(); } @@ -115,12 +107,12 @@ template void static_matrix::init_empty_matrix init_row_columns(m, n); } -template unsigned static_matrix::lowest_row_in_column(unsigned col) { +template unsigned static_matrix::lowest_row_in_column(unsigned col) { + lp_assert(col < column_count()); column_strip & colstrip = m_columns[col]; - lp_assert(colstrip.live_size() > 0); + lp_assert(colstrip.size() > 0); unsigned ret = 0; - for (const auto & t : colstrip) { - if (t.dead()) continue; + for (auto & t : colstrip) { if (t.var() > ret) { ret = t.var(); } @@ -147,7 +139,7 @@ template void static_matrix::remove_last_column(u 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->var() == j) { + if (row_it-.var() == j) { row.erase(row.begin() + offset); break; } @@ -158,6 +150,18 @@ 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; @@ -262,9 +266,50 @@ 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].var() > k) { + col[i].var()--; + } + } + } +} + +template void static_matrix::cross_out_row_from_columns(unsigned k, row_strip & row) { + for (auto & t : row) { + cross_out_row_from_column(t.var(), 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].var() == 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.dead()) continue; if (t.var() == j) { return t.get_val(); } @@ -297,22 +342,16 @@ template bool static_matrix::is_correct() const { auto &r = m_rows[i]; std::unordered_set s; for (auto & rc : r) { - if (rc.dead()) continue; if (s.find(rc.var()) != s.end()) { return false; } s.insert(rc.var()); if (rc.var() >= m_columns.size()) return false; - const auto& col = m_columns[rc.var()]; - if (col.cells_size() <= rc.offset()) + if (rc.m_offset >= m_columns[rc.var()].size()) return false; - const auto &cc = col[rc.offset()]; - if (cc.dead()) + if (rc.get_val() != get_val(m_columns[rc.var()][rc.m_offset])) return false; - if (& m_rows[cc.var()][cc.offset()] != & rc) { - return false; - } if (is_zero(rc.get_val())) { return false; } @@ -324,56 +363,51 @@ template bool static_matrix::is_correct() const { auto & c = m_columns[j]; std::unordered_set s; for (auto & cc : c) { - if (cc.dead()) - continue; if (s.find(cc.var()) != s.end()) { return false; } s.insert(cc.var()); if (cc.var() >= m_rows.size()) return false; - auto & rc = m_rows[cc.var()][cc.offset()]; - if (rc.dead()) + if (cc.m_offset >= m_rows[cc.var()].size()) return false; - if (&cc != &m_columns[rc.var()][rc.offset()]) + if (get_val(cc) != m_rows[cc.var()][cc.m_offset].get_val()) 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(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); -} +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.var()]; + column_cell& cs = m_columns[row_el_iv.var()][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.var()][cc.offset()].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.var()][rc.offset()].offset() = row_offset; + } + column_vals.pop_back(); + row_vals.pop_back(); +} template -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; +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)); } } From 7370396c3052d06cb8e47f4eb062ea6cca9bd53b Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Wed, 1 Aug 2018 08:06:56 -0700 Subject: [PATCH 2/4] fix the build Signed-off-by: Lev Nachmanson --- src/util/lp/static_matrix_def.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/lp/static_matrix_def.h b/src/util/lp/static_matrix_def.h index 1dcaa34f8..7949573de 100644 --- a/src/util/lp/static_matrix_def.h +++ b/src/util/lp/static_matrix_def.h @@ -139,7 +139,7 @@ template void static_matrix::remove_last_column(u 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-.var() == j) { + if (row_it.var() == j) { row.erase(row.begin() + offset); break; } From 075cf106aae7285df7b89cfdfcf3549d487648bf Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Wed, 1 Aug 2018 08:46:03 -0700 Subject: [PATCH 3/4] fix the build Signed-off-by: Lev Nachmanson --- src/util/lp/lar_solver.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/util/lp/lar_solver.cpp b/src/util/lp/lar_solver.cpp index 78bd45e94..654eb7017 100644 --- a/src/util/lp/lar_solver.cpp +++ b/src/util/lp/lar_solver.cpp @@ -870,10 +870,10 @@ void lar_solver::fill_last_row_of_A_r(static_matrix> & A, for (auto & t : ls->m_coeffs) { lp_assert(!is_zero(t.second)); var_index j = t.first; - A.add_new_element(last_row, j, - t.second); + A.set(last_row, j, - t.second); } unsigned basis_j = A.column_count() - 1; - A.add_new_element(last_row, basis_j, mpq(1)); + A.set(last_row, basis_j, mpq(1)); } template @@ -895,7 +895,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.add_new_element(i, it.var(), convert_struct::convert(it.get_val())); + matr.set(i, it.var(), convert_struct::convert(it.get_val())); } } } @@ -1840,11 +1840,11 @@ 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.add_new_element(last_row, j, -t.second.get_double()); + A.set(last_row, j, -t.second.get_double()); } unsigned basis_j = A.column_count() - 1; - A.add_new_element(last_row, basis_j, -1); + A.set(last_row, basis_j, -1); lp_assert(A.is_correct()); } From adfbc2d0012b7ff8538eea2599ee3a05be28b4b7 Mon Sep 17 00:00:00 2001 From: Lev Nachmanson Date: Wed, 1 Aug 2018 10:26:39 -0700 Subject: [PATCH 4/4] fix the build Signed-off-by: Lev Nachmanson --- src/test/lp/lp.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp index 37a096827..6e418fe68 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.add_new_element(i, j, v); + A.set(i, j, v); } while (true); }