mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-15 13:28:47 +00:00
Code Activity

work on static_matrix's cells

Browse source 
Signed-off-by: Lev <levnach@hotmail.com>

trying the new scheme in static_matrix : in progress

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

in the middle of changes in static_matrix

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

more fixes in static_matrix.h

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

debug

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fixes in static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

fixes in static_matrix, column_strip

Signed-off-by: Lev <levnach@hotmail.com>

fixes in static_matrix

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

fixes for static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

work on static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

work on static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

progress in static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

fix a bug in swap_with_head_cell

Signed-off-by: Lev <levnach@hotmail.com>

progress in static_matrix

Signed-off-by: Lev <levnach@hotmail.com>

compress rows and columns if needed

Signed-off-by: Lev <levnach@hotmail.com>

fix in compression of cells

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2018-07-18 12:50:54 -07:00 committed by Lev
parent 95845bbb01
commit 16b71fe911
22 changed files with 767 additions and 388 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
src/shell/main.cpp
View file

@ -295,7 +295,13 @@ void parse_cmd_line_args(int argc, char ** argv) {
int STD_CALL main(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; unsigned return_value = 0;
memory::initialize(0); memory::initialize(0);
memory::exit_when_out_of_memory(true, "ERROR: out of memory"); memory::exit_when_out_of_memory(true, "ERROR: out of memory");

4
src/tactic/smtlogics/qflia_tactic.cpp
View file

@ -212,6 +212,9 @@ tactic * mk_qflia_tactic(ast_manager & m, params_ref const & p) {
tactic * st = using_params(and_then(preamble_st, tactic * st = using_params(and_then(preamble_st,
#if 1
mk_smt_tactic()),
#else
or_else(mk_ilp_model_finder_tactic(m), or_else(mk_ilp_model_finder_tactic(m),
mk_pb_tactic(m), mk_pb_tactic(m),
and_then(fail_if_not(mk_is_quasi_pb_probe()), 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_fail_if_undecided_tactic()),
mk_bounded_tactic(m), mk_bounded_tactic(m),
mk_smt_tactic())), mk_smt_tactic())),
#endif
main_p); main_p);

2
src/test/lp/lp.cpp
View file

@ -2541,7 +2541,7 @@ void read_row_cols(unsigned i, static_matrix<double, double>& A, std::ifstream &
lp_assert(r.size() == 4); lp_assert(r.size() == 4);
unsigned j = atoi(r[1].c_str()); unsigned j = atoi(r[1].c_str());
double v = atof(r[3].c_str()); double v = atof(r[3].c_str());
A.set(i, j, v); A.add_new_element(i, j, v);
} while (true); } while (true);
} }

70
src/util/lp/bound_analyzer_on_row.h
View file

@ -29,65 +29,7 @@ Revision History:
namespace lp { namespace lp {
template <typename C> // C plays a role of a container template <typename C> // C plays a role of a container
class bound_analyzer_on_row { class bound_analyzer_on_row {
struct term_with_basis_col { const C& m_row;
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<unsigned>(-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<unsigned>(-1))
m_it++;
else
m_bj = static_cast<unsigned>(-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;
bound_propagator & m_bp; bound_propagator & m_bp;
unsigned m_row_or_term_index; unsigned m_row_or_term_index;
int m_column_of_u; // index of an unlimited from above monoid int m_column_of_u; // index of an unlimited from above monoid
@ -105,7 +47,7 @@ public :
bound_propagator & bp bound_propagator & bp
) )
: :
m_row(it, bj), m_row(it),
m_bp(bp), m_bp(bp),
m_row_or_term_index(row_or_term_index), m_row_or_term_index(row_or_term_index),
m_column_of_u(-1), m_column_of_u(-1),
@ -117,6 +59,8 @@ public :
unsigned j; unsigned j;
void analyze() { void analyze() {
for (const auto & c : m_row) { for (const auto & c : m_row) {
if (c.dead())
continue;
if ((m_column_of_l == -2) && (m_column_of_u == -2)) if ((m_column_of_l == -2) && (m_column_of_u == -2))
break; break;
analyze_bound_on_var_on_coeff(c.var(), c.coeff()); analyze_bound_on_var_on_coeff(c.var(), c.coeff());
@ -226,6 +170,7 @@ public :
mpq total; mpq total;
lp_assert(is_zero(total)); lp_assert(is_zero(total));
for (const auto& p : m_row) { for (const auto& p : m_row) {
if (p.dead()) continue;
bool str; bool str;
total -= monoid_min(p.coeff(), p.var(), str); total -= monoid_min(p.coeff(), p.var(), str);
if (str) if (str)
@ -234,6 +179,7 @@ public :
for (const auto &p : m_row) { for (const auto &p : m_row) {
if (p.dead()) continue;
bool str; bool str;
bool a_is_pos = is_pos(p.coeff()); bool a_is_pos = is_pos(p.coeff());
mpq bound = total / p.coeff() + monoid_min_no_mult(a_is_pos, p.var(), str); mpq bound = total / p.coeff() + monoid_min_no_mult(a_is_pos, p.var(), str);
@ -251,6 +197,7 @@ public :
mpq total; mpq total;
lp_assert(is_zero(total)); lp_assert(is_zero(total));
for (const auto &p : m_row) { for (const auto &p : m_row) {
if (p.dead()) continue;
bool str; bool str;
total -= monoid_max(p.coeff(), p.var(), str); total -= monoid_max(p.coeff(), p.var(), str);
if (str) if (str)
@ -258,6 +205,7 @@ public :
} }
for (const auto& p : m_row) { for (const auto& p : m_row) {
if (p.dead()) continue;
bool str; bool str;
bool a_is_pos = is_pos(p.coeff()); bool a_is_pos = is_pos(p.coeff());
mpq bound = total / p.coeff() + monoid_max_no_mult(a_is_pos, p.var(), str); 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; mpq bound = -m_rs.x;
bool strict = false; bool strict = false;
for (const auto& p : m_row) { for (const auto& p : m_row) {
if (p.dead()) continue;
j = p.var(); j = p.var();
if (j == static_cast<unsigned>(m_column_of_u)) { if (j == static_cast<unsigned>(m_column_of_u)) {
u_coeff = p.coeff(); u_coeff = p.coeff();
@ -309,6 +258,7 @@ public :
mpq bound = -m_rs.x; mpq bound = -m_rs.x;
bool strict = false; bool strict = false;
for (const auto &p : m_row) { for (const auto &p : m_row) {
if (p.dead()) continue;
j = p.var(); j = p.var();
if (j == static_cast<unsigned>(m_column_of_l)) { if (j == static_cast<unsigned>(m_column_of_l)) {
l_coeff = p.coeff(); l_coeff = p.coeff();

2
src/util/lp/bound_propagator.cpp
View file

@ -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]; 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) { 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); j = m_lar_solver.adjust_column_index_to_term_index(j);
if (m_lar_solver.is_term(j)) { if (m_lar_solver.is_term(j)) {
// lp treats terms as not having a free coefficient, restoring it below for the outside consumption // lp treats terms as not having a free coefficient, restoring it below for the outside consumption

2
src/util/lp/core_solver_pretty_printer_def.h
View file

@ -100,7 +100,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_
for (unsigned column = 0; column < ncols(); column++) { for (unsigned column = 0; column < ncols(); column++) {
vector<T> t(nrows(), zero_of_type<T>()); vector<T> t(nrows(), zero_of_type<T>());
for (const auto & c : m_core_solver.m_A.m_columns[column]){ 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); string name = m_core_solver.column_name(column);

1
src/util/lp/indexed_vector.h
View file

@ -172,6 +172,7 @@ public:
} }
unsigned var() const { return m_var;} unsigned var() const { return m_var;}
const T & coeff() const { return m_coeff; } const T & coeff() const { return m_coeff; }
bool dead() const { return false; }
}; };
struct const_iterator { struct const_iterator {

17
src/util/lp/int_solver.cpp
View file

@ -100,6 +100,7 @@ bool int_solver::is_gomory_cut_target(const row_strip<mpq>& row) {
// All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed). // All non base variables must be at their bounds and assigned to rationals (that is, infinitesimals are not allowed).
unsigned j; unsigned j;
for (const auto & p : row) { for (const auto & p : row) {
if (p.dead()) continue;
j = p.var(); j = p.var();
if (is_base(j)) continue; if (is_base(j)) continue;
if (!at_bound(j)) if (!at_bound(j))
@ -311,6 +312,7 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip<mpq> & row
TRACE("gomory_cut", TRACE("gomory_cut",
tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl; tout << "applying cut at:\n"; m_lar_solver->print_row(row, tout); tout << std::endl;
for (auto & p : row) { 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); m_lar_solver->m_mpq_lar_core_solver.m_r_solver.print_column_info(p.var(), tout);
} }
tout << "inf_col = " << inf_col << std::endl; tout << "inf_col = " << inf_col << std::endl;
@ -324,7 +326,8 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip<mpq> & row
bool some_int_columns = false; bool some_int_columns = false;
mpq f_0 = int_solver::fractional_part(get_value(inf_col)); mpq f_0 = int_solver::fractional_part(get_value(inf_col));
mpq one_min_f_0 = 1 - f_0; 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(); x_j = p.var();
if (x_j == inf_col) if (x_j == inf_col)
continue; continue;
@ -353,6 +356,7 @@ lia_move int_solver::mk_gomory_cut( unsigned inf_col, const row_strip<mpq> & row
int int_solver::find_free_var_in_gomory_row(const row_strip<mpq>& row) { int int_solver::find_free_var_in_gomory_row(const row_strip<mpq>& row) {
unsigned j; unsigned j;
for (const auto & p : row) { for (const auto & p : row) {
if (p.dead()) continue;
j = p.var(); j = p.var();
if (!is_base(j) && is_free(j)) if (!is_base(j) && is_free(j))
return static_cast<int>(j); return static_cast<int>(j);
@ -789,6 +793,7 @@ lia_move int_solver::patch_nbasic_columns() {
mpq get_denominators_lcm(const row_strip<mpq> & row) { mpq get_denominators_lcm(const row_strip<mpq> & row) {
mpq r(1); mpq r(1);
for (auto & c : row) { for (auto & c : row) {
if (c.dead()) continue;
r = lcm(r, denominator(c.coeff())); r = lcm(r, denominator(c.coeff()));
} }
return r; return r;
@ -802,6 +807,7 @@ bool int_solver::gcd_test_for_row(static_matrix<mpq, numeric_pair<mpq>> & A, uns
bool least_coeff_is_bounded = false; bool least_coeff_is_bounded = false;
unsigned j; unsigned j;
for (auto &c : A.m_rows[i]) { for (auto &c : A.m_rows[i]) {
if (c.dead()) continue;
j = c.var(); j = c.var();
const mpq& a = c.coeff(); const mpq& a = c.coeff();
if (m_lar_solver->column_is_fixed(j)) { 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<mpq> & row) { void int_solver::fill_explanation_from_fixed_columns(const row_strip<mpq> & row) {
for (const auto & c : row) { for (const auto & c : row) {
if (c.dead()) continue;
if (!m_lar_solver->column_is_fixed(c.var())) if (!m_lar_solver->column_is_fixed(c.var()))
continue; continue;
add_to_explanation_from_fixed_or_boxed_column(c.var()); add_to_explanation_from_fixed_or_boxed_column(c.var());
@ -892,6 +899,7 @@ bool int_solver::ext_gcd_test(const row_strip<mpq> & row,
mpq a; mpq a;
unsigned j; unsigned j;
for (const auto & c : row) { for (const auto & c : row) {
if (c.dead()) continue;
j = c.var(); j = c.var();
const mpq & a = c.coeff(); const mpq & a = c.coeff();
if (m_lar_solver->column_is_fixed(j)) 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()); lp_assert(settings().use_tableau());
const auto & A = m_lar_solver->A_r(); const auto & A = m_lar_solver->A_r();
for (const auto &c : A.column(j)) { for (const auto &c : A.column(j)) {
if (c.dead()) continue;
row_index = c.var(); row_index = c.var();
const mpq & a = c.coeff(); 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 { 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; 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<mpq>::is_pos(c.coeff())) if (numeric_traits<mpq>::is_pos(c.coeff()))
out << "+"; out << "+";
out << c.coeff() << rslv.column_name(c.var()) << " "; 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(c.var(), out);
} }
rslv.print_column_bound_info(rslv.m_basis[row_index], out); rslv.print_column_bound_info(rslv.m_basis[row_index], out);

10
src/util/lp/lar_core_solver.h
View file

@ -339,7 +339,7 @@ public:
if (!update_xj_and_get_delta(j, pos_type, delta)) if (!update_xj_and_get_delta(j, pos_type, delta))
continue; continue;
for (const auto & cc : m_r_solver.m_A.m_columns[j]){ 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]; 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)); 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()) if (!m_r_solver.m_settings.use_tableau())
return true; return true;
for (unsigned j : m_r_solver.m_basis) { 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.m_columns[j].live_size() == 1);
lp_assert(m_r_solver.m_A.get_val(m_r_solver.m_A.m_columns[j][0]) == one_of_type<mpq>());
} }
for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) { for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) {
if (m_r_solver.m_basis_heading[j] >= 0) continue; if (m_r_solver.m_basis_heading[j] >= 0) continue;
@ -633,8 +632,9 @@ public:
void create_double_matrix(static_matrix<double, double> & A) { void create_double_matrix(static_matrix<double, double> & A) {
for (unsigned i = 0; i < m_r_A.row_count(); i++) { for (unsigned i = 0; i < m_r_A.row_count(); i++) {
auto & row = m_r_A.m_rows[i]; auto & row = m_r_A.m_rows[i];
for (row_cell<mpq> & c : row) { for (row_cell<mpq> & c : row.m_cells) {
A.set(i, c.m_j, c.get_val().get_double()); if (c.dead()) continue;
A.add_new_element(i, c.var(), c.get_val().get_double());
} }
} }
} }

5
src/util/lp/lar_core_solver_def.h
View file

@ -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]; 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_sum_sign = m_r_solver.inf_sign_of_column(bj);
m_infeasible_linear_combination.clear(); m_infeasible_linear_combination.clear();
for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) { for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau].m_cells) {
m_infeasible_linear_combination.push_back(std::make_pair( rc.get_val(), rc.m_j)); if (rc.dead()) continue;
m_infeasible_linear_combination.push_back(std::make_pair( rc.get_val(), rc.var()));
} }
} }

66
src/util/lp/lar_solver.cpp
View file

@ -162,9 +162,8 @@ void lar_solver::analyze_new_bounds_on_row_tableau(
unsigned row_index, unsigned row_index,
bound_propagator & bp ) { 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; return;
lp_assert(use_tableau()); lp_assert(use_tableau());
bound_analyzer_on_row<row_strip<mpq>>::analyze_row(A_r().m_rows[row_index], bound_analyzer_on_row<row_strip<mpq>>::analyze_row(A_r().m_rows[row_index],
static_cast<unsigned>(-1), static_cast<unsigned>(-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); bound_j = m_var_register.external_to_local(bound_j);
} }
for (auto const& r : A_r().m_rows[i]) { 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; if (j == bound_j) continue;
mpq const& a = r.get_val(); mpq const& a = r.get_val();
int a_sign = is_pos(a)? 1: -1; 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) if (i < 0)
jset.insert(j); jset.insert(j);
else { else {
for (auto & rc : A_r().m_rows[i]) for (const auto & rc : A_r().m_rows[i])
jset.insert(rc.m_j); 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) { 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]) 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(); } 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 { bool lar_solver::row_is_correct(unsigned i) const {
numeric_pair<mpq> r = zero_of_type<numeric_pair<mpq>>(); numeric_pair<mpq> r = zero_of_type<numeric_pair<mpq>>();
for (const auto & c : A_r().m_rows[i]) for (const auto & c : A_r().m_rows[i]) {
r += c.m_value * m_mpq_lar_core_solver.m_r_x[c.m_j]; if (c.dead()) continue;
r += c.coeff() * m_mpq_lar_core_solver.m_r_x[c.var()];
}
return is_zero(r); 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<mpq> & delta) { void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq> & delta) {
if (use_tableau()) { if (use_tableau()) {
for (const auto & c : A_r().m_columns[j]) { 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()) { if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj); m_basic_columns_with_changed_cost.insert(bj);
} }
@ -790,10 +795,10 @@ numeric_pair<mpq> lar_solver::get_basic_var_value_from_row_directly(unsigned i)
unsigned bj = m_mpq_lar_core_solver.m_r_solver.m_basis[i]; unsigned bj = m_mpq_lar_core_solver.m_r_solver.m_basis[i];
for (const auto & c: A_r().m_rows[i]) { for (const auto & c: A_r().m_rows[i]) {
if (c.m_j == bj) continue; if (c.var() == bj) continue;
const auto & x = m_mpq_lar_core_solver.m_r_x[c.m_j]; const auto & x = m_mpq_lar_core_solver.m_r_x[c.var()];
if (!is_zero(x)) if (!is_zero(x))
r -= c.m_value * x; r -= c.coeff() * x;
} }
return r; return r;
} }
@ -866,14 +871,14 @@ void lar_solver::fill_last_row_of_A_r(static_matrix<mpq, numeric_pair<mpq>> & A,
lp_assert(A.row_count() > 0); lp_assert(A.row_count() > 0);
lp_assert(A.column_count() > 0); lp_assert(A.column_count() > 0);
unsigned last_row = A.row_count() - 1; 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) { for (auto & t : ls->m_coeffs) {
lp_assert(!is_zero(t.second)); lp_assert(!is_zero(t.second));
var_index j = t.first; 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; 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 <typename U, typename V> template <typename U, typename V>
@ -895,7 +900,7 @@ void lar_solver::copy_from_mpq_matrix(static_matrix<U, V> & matr) {
matr.m_columns.resize(A_r().column_count()); matr.m_columns.resize(A_r().column_count());
for (unsigned i = 0; i < matr.row_count(); i++) { for (unsigned i = 0; i < matr.row_count(); i++) {
for (auto & it : A_r().m_rows[i]) { for (auto & it : A_r().m_rows[i]) {
matr.set(i, it.m_j, convert_struct<U, mpq>::convert(it.get_val())); matr.add_new_element(i, it.var(), convert_struct<U, mpq>::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); lp_assert(A_r().row_count() == i + 1 && A_r().column_count() == j + 1);
auto & last_column = A_r().m_columns[j]; auto & last_column = A_r().m_columns[j];
int non_zero_column_cell_index = -1; 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]; auto & cc = last_column[k];
if (cc.m_i == i) if (cc.dead()) continue;
if (cc.var() == i)
return; return;
non_zero_column_cell_index = k; non_zero_column_cell_index = k;
} }
lp_assert(non_zero_column_cell_index != -1); lp_assert(non_zero_column_cell_index != -1);
lp_assert(static_cast<unsigned>(non_zero_column_cell_index) != i); lp_assert(static_cast<unsigned>(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) { 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]; auto & last_row = A_r().m_rows[i];
mpq &cost_j = m_mpq_lar_core_solver.m_r_solver.m_costs[j]; mpq &cost_j = m_mpq_lar_core_solver.m_r_solver.m_costs[j];
bool cost_is_nz = !is_zero(cost_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]; auto &rc = last_row[k];
if (rc.dead()) {
last_row.pop();
continue;
}
if (cost_is_nz) { 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(last_row.live_size() == 0);
lp_assert(A_r().m_columns[j].size() == 0); lp_assert(A_r().m_columns[j].live_size() == 0);
A_r().m_rows.pop_back(); A_r().m_rows.pop_back();
A_r().m_columns.pop_back(); A_r().m_columns.pop_back();
lp_assert(A_r().is_correct());
slv.m_b.pop_back(); slv.m_b.pop_back();
} }
void lar_solver::remove_last_column_from_A() { void lar_solver::remove_last_column_from_A() {
// the last column has to be empty // 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(); A_r().m_columns.pop_back();
} }
@ -1840,11 +1851,12 @@ void lar_solver::fill_last_row_of_A_d(static_matrix<double, double> & A, const l
for (auto & t : ls->m_coeffs) { for (auto & t : ls->m_coeffs) {
lp_assert(!is_zero(t.second)); lp_assert(!is_zero(t.second));
var_index j = t.first; 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; 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) { void lar_solver::update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {

18
src/util/lp/lp_core_solver_base.h
View file

@ -213,8 +213,17 @@ public:
unsigned m = m_A.row_count(); unsigned m = m_A.row_count();
for (unsigned i = 0; i < m; i++) { for (unsigned i = 0; i < m; i++) {
unsigned bj = m_basis[i]; unsigned bj = m_basis[i];
lp_assert(m_A.m_columns[bj].size() > 0); lp_assert(m_A.m_columns[bj].live_size() > 0);
if (m_A.m_columns[bj].size() > 1 || m_A.get_val(m_A.m_columns[bj][0]) != one_of_type<mpq>()) return true; 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<mpq>())
return true;
else
break;
}
} }
return false; return false;
} }
@ -237,8 +246,9 @@ public:
} }
} else { } else {
auto d = m_costs[j]; auto d = m_costs[j];
for (auto & cc : this->m_A.m_columns[j]) { for (const auto & cc : this->m_A.m_columns[j]) {
d -= this->m_costs[this->m_basis[cc.m_i]] * this->m_A.get_val(cc); if (cc.dead()) continue;
d -= this->m_costs[this->m_basis[cc.var()]] * this->m_A.get_val(cc);
} }
if (m_d[j] != d) { if (m_d[j] != d) {
return false; return false;

91
src/util/lp/lp_core_solver_base_def.h
View file

@ -103,9 +103,12 @@ pivot_to_reduced_costs_tableau(unsigned i, unsigned j) {
T &a = m_d[j]; T &a = m_d[j];
if (is_zero(a)) if (is_zero(a))
return; return;
for (const row_cell<T> & r: m_A.m_rows[i]) for (const row_cell<T> & r: m_A.m_rows[i]){
if (r.m_j != j) if (r.dead()) continue;
m_d[r.m_j] -= a * r.get_val(); if (r.var() != j)
m_d[r.var()] -= a * r.get_val();
}
ls
a = zero_of_type<T>(); // zero the pivot column's m_d finally a = zero_of_type<T>(); // 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<T>::is_zero(pi_1)) { if (numeric_traits<T>::is_zero(pi_1)) {
continue; continue;
} }
for (auto & c : m_A.m_rows[i]) { for (auto & c : m_A.m_rows[i].m_cells) {
unsigned j = c.m_j; if (c.dead()) continue;
unsigned j = c.var();
if (m_basis_heading[j] < 0) { if (m_basis_heading[j] < 0) {
m_pivot_row.add_value_at_index_with_drop_tolerance(j, c.get_val() * pi_1); 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 else
for (const auto & c : m_A.m_columns[entering]) { 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); m_x[m_basis[i]] -= delta * m_A.get_val(c);
} }
} }
@ -449,7 +453,7 @@ rs_minus_Anx(vector<X> & rs) {
while (row--) { while (row--) {
auto &rsv = rs[row] = m_b[row]; auto &rsv = rs[row] = m_b[row];
for (auto & it : m_A.m_rows[row]) { for (auto & it : m_A.m_rows[row]) {
unsigned j = it.m_j; unsigned j = it.var();
if (m_basis_heading[j] < 0) { if (m_basis_heading[j] < 0) {
rsv -= m_x[j] * it.get_val(); 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<T>::precise()); lp_assert(numeric_traits<T>::precise());
int pivot_index = -1; int pivot_index = -1;
auto & row = m_A.m_rows[pivot_row]; auto & row = m_A.m_rows[pivot_row];
unsigned size = row.size(); unsigned size = row.cells_size();
for (unsigned j = 0; j < size; j++) { 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<int>(j); pivot_index = static_cast<int>(j);
break; break;
} }
@ -599,26 +605,33 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) {
if (pivot_index == -1) if (pivot_index == -1)
return false; return false;
auto & pivot_cell = row[pivot_index]; auto & pivot_cell = row[pivot_index];
if (is_zero(pivot_cell.m_value)) T & coeff = pivot_cell.coeff();
if (is_zero(coeff))
return false; return false;
this->m_b[pivot_row] /= pivot_cell.m_value; this->m_b[pivot_row] /= coeff;
for (unsigned j = 0; j < size; j++) { for (unsigned j = 0; j < size; j++) {
if (row[j].m_j != pivot_col) { auto & c = row[j];
row[j].m_value /= pivot_cell.m_value; if (c.dead()) continue;
if (c.var() != pivot_col) {
c.coeff() /= coeff;
} }
} }
pivot_cell.m_value = one_of_type<T>(); coeff = one_of_type<T>();
lp_assert(m_A.is_correct());
return true; return true;
} }
template <typename T, typename X> bool lp_core_solver_base<T, X>:: template <typename T, typename X> bool lp_core_solver_base<T, X>::
pivot_column_tableau(unsigned j, unsigned piv_row_index) { 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; return false;
auto &column = m_A.m_columns[j]; auto &column = m_A.m_columns[j];
int pivot_col_cell_index = -1; int pivot_col_cell_index = -1;
for (unsigned k = 0; k < column.size(); k++) { for (unsigned k = 0; k < column.cells_size(); k++) {
if (column[k].m_i == piv_row_index) { if (column[k].dead()) continue;
if (column[k].var() == piv_row_index) {
pivot_col_cell_index = k; pivot_col_cell_index = k;
break; break;
} }
@ -626,25 +639,26 @@ pivot_column_tableau(unsigned j, unsigned piv_row_index) {
if (pivot_col_cell_index < 0) if (pivot_col_cell_index < 0)
return false; return false;
if (pivot_col_cell_index != 0) { if (pivot_col_cell_index != 0)
lp_assert(column.size() > 1); m_A.swap_with_head_cell(j, pivot_col_cell_index);
// 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;
m_A.m_rows[piv_row_index][column[0].m_offset].m_offset = 0; lp_assert(m_A.is_correct());
m_A.m_rows[c.m_i][c.m_offset].m_offset = pivot_col_cell_index; while (column.live_size() > 1) {
}
while (column.size() > 1) {
auto & c = column.back(); 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)) { if(! m_A.pivot_row_to_row_given_cell(piv_row_index, c, j)) {
return false; return false;
} }
if (m_pivoted_rows!= nullptr) 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) if (m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs)
pivot_to_reduced_costs_tableau(piv_row_index, j); pivot_to_reduced_costs_tableau(piv_row_index, j);
@ -758,10 +772,11 @@ fill_reduced_costs_from_m_y_by_rows() {
while (i--) { while (i--) {
const T & y = m_y[i]; const T & y = m_y[i];
if (is_zero(y)) continue; if (is_zero(y)) continue;
for (row_cell<T> & it : m_A.m_rows[i]) { for (row_cell<T> & c : m_A.m_rows[i].m_cells) {
j = it.m_j; if (c.dead()) continue;
j = c.var();
if (m_basis_heading[j] < 0) { 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<X>& rs){
while (row--) { while (row--) {
auto &rsv = rs[row]; auto &rsv = rs[row];
for (auto & it : m_A.m_rows[row]) { for (auto & it : m_A.m_rows[row]) {
unsigned j = it.m_j; unsigned j = it.var();
if (m_basis_heading[j] >= 0) { if (m_basis_heading[j] >= 0) {
rsv -= m_x[j] * it.get_val(); rsv -= m_x[j] * it.get_val();
} }
@ -957,7 +972,8 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::pivot_fixed_v
if (get_column_type(basic_j) != column_type::fixed) continue; if (get_column_type(basic_j) != column_type::fixed) continue;
T a; T a;
unsigned j; 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(); j = c.var();
if (j == basic_j) if (j == basic_j)
continue; continue;
@ -972,7 +988,8 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::pivot_fixed_v
template <typename T, typename X> bool lp_core_solver_base<T, X>::remove_from_basis(unsigned basic_j) { template <typename T, typename X> bool lp_core_solver_base<T, X>::remove_from_basis(unsigned basic_j) {
indexed_vector<T> w(m_basis.size()); // the buffer indexed_vector<T> w(m_basis.size()); // the buffer
unsigned i = m_basis_heading[basic_j]; 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) if (c.var() == basic_j)
continue; continue;
if (pivot_column_general(c.var(), basic_j, w)) if (pivot_column_general(c.var(), basic_j, w))
@ -1043,8 +1060,8 @@ void lp_core_solver_base<T, X>::calculate_pivot_row(unsigned i) {
if (m_settings.use_tableau()) { if (m_settings.use_tableau()) {
unsigned basic_j = m_basis[i]; unsigned basic_j = m_basis[i];
for (auto & c : m_A.m_rows[i]) { for (auto & c : m_A.m_rows[i]) {
if (c.m_j != basic_j) if (c.var() != basic_j)
m_pivot_row.set_value(c.get_val(), c.m_j); m_pivot_row.set_value(c.get_val(), c.var());
} }
return; return;
} }

62
src/util/lp/lp_primal_core_solver.h
View file

@ -172,7 +172,7 @@ public:
bool monoid_can_decrease(const row_cell<T> & rc) const { bool monoid_can_decrease(const row_cell<T> & rc) const {
unsigned j = rc.m_j; unsigned j = rc.var();
lp_assert(this->column_is_feasible(j)); lp_assert(this->column_is_feasible(j));
switch (this->m_column_types[j]) { switch (this->m_column_types[j]) {
case column_type::free_column: case column_type::free_column:
@ -205,7 +205,7 @@ public:
} }
bool monoid_can_increase(const row_cell<T> & rc) const { bool monoid_can_increase(const row_cell<T> & rc) const {
unsigned j = rc.m_j; unsigned j = rc.var();
lp_assert(this->column_is_feasible(j)); lp_assert(this->column_is_feasible(j));
switch (this->m_column_types[j]) { switch (this->m_column_types[j]) {
case column_type::free_column: 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 get_number_of_basic_vars_that_might_become_inf(unsigned j) const { // consider looking at the signs here: todo
unsigned r = 0; unsigned r = 0;
for (auto & cc : this->m_A.m_columns[j]) { for (const auto & cc : this->m_A.m_columns[j]) {
unsigned k = this->m_basis[cc.m_i]; if (cc.dead()) continue;
unsigned k = this->m_basis[cc.var()];
if (this->m_column_types[k] != column_type::free_column) if (this->m_column_types[k] != column_type::free_column)
r++; r++;
} }
@ -253,7 +254,8 @@ public:
unsigned bj = this->m_basis[i]; unsigned bj = this->m_basis[i];
bool bj_needs_to_grow = needs_to_grow(bj); bool bj_needs_to_grow = needs_to_grow(bj);
for (const row_cell<T>& rc : this->m_A.m_rows[i]) { for (const row_cell<T>& rc : this->m_A.m_rows[i]) {
if (rc.m_j == bj) if (rc.dead()) continue;
if (rc.var() == bj)
continue; continue;
if (bj_needs_to_grow) { if (bj_needs_to_grow) {
if (!monoid_can_decrease(rc)) if (!monoid_can_decrease(rc))
@ -262,9 +264,9 @@ public:
if (!monoid_can_increase(rc)) if (!monoid_can_increase(rc))
continue; continue;
} }
if (rc.m_j < static_cast<unsigned>(j) ) { if (rc.var() < static_cast<unsigned>(j) ) {
j = rc.m_j; j = rc.var();
a_ent = rc.m_value; a_ent = rc.coeff();
} }
} }
if (j == -1) { if (j == -1) {
@ -278,13 +280,18 @@ public:
if (m_bland_mode_tableau) if (m_bland_mode_tableau)
return find_beneficial_column_in_row_tableau_rows_bland_mode(i, a_ent); 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 // a short row produces short infeasibility explanation and benefits at least one pivot operation
vector<const row_cell<T>*> choices; int choice = -1;
unsigned num_of_non_free_basics = 1000000; int nchoices = 0;
unsigned len = 100000000; unsigned len = 100000000;
unsigned bj = this->m_basis[i]; unsigned bj = this->m_basis[i];
bool bj_needs_to_grow = needs_to_grow(bj); bool bj_needs_to_grow = needs_to_grow(bj);
for (const row_cell<T>& rc : this->m_A.m_rows[i]) { for (unsigned k = 0; k < this->m_A.m_rows[i].m_cells.size(); k++) {
unsigned j = rc.m_j; const row_cell<T>& rc = this->m_A.m_rows[i].m_cells[k];
if (rc.dead()) continue;
const row_cell<T>& 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) if (j == bj)
continue; continue;
if (bj_needs_to_grow) { if (bj_needs_to_grow) {
@ -297,26 +304,24 @@ public:
unsigned damage = get_number_of_basic_vars_that_might_become_inf(j); unsigned damage = get_number_of_basic_vars_that_might_become_inf(j);
if (damage < num_of_non_free_basics) { if (damage < num_of_non_free_basics) {
num_of_non_free_basics = damage; num_of_non_free_basics = damage;
len = this->m_A.m_columns[j].size(); len = this->m_A.m_columns[j].live_size();
choices.clear(); choice = k;
choices.push_back(&rc); nchoices = 1;
} else if (damage == num_of_non_free_basics && } else if (damage == num_of_non_free_basics &&
this->m_A.m_columns[j].size() <= len && (this->m_settings.random_next() % 2)) { this->m_A.m_columns[j].live_size() <= len && (this->m_settings.random_next() % (++nchoices))) {
choices.push_back(&rc); choice = k;
len = this->m_A.m_columns[j].size(); len = this->m_A.m_columns[j].live_size();
} }
} }
if (choices.size() == 0) { if (choice == -1) {
m_inf_row_index_for_tableau = i; m_inf_row_index_for_tableau = i;
return -1; return -1;
} }
const row_cell<T>* rc = choices.size() == 1? choices[0] : const row_cell<T>& rc = this->m_A.m_rows[i].m_cells[choice];
choices[this->m_settings.random_next() % choices.size()]; a_ent = rc.coeff();
return rc.var();
a_ent = rc->m_value;
return rc->m_j;
} }
static X positive_infinity() { static X positive_infinity() {
return convert_struct<X, unsigned>::convert(std::numeric_limits<unsigned>::max()); return convert_struct<X, unsigned>::convert(std::numeric_limits<unsigned>::max());
@ -827,11 +832,11 @@ public:
this->m_rows_nz.resize(this->m_A.row_count()); this->m_rows_nz.resize(this->m_A.row_count());
for (unsigned i = 0; i < this->m_A.column_count(); i++) { for (unsigned i = 0; i < this->m_A.column_count(); i++) {
if (this->m_columns_nz[i] == 0) 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++) { for (unsigned i = 0; i < this->m_A.row_count(); i++) {
if (this->m_rows_nz[i] == 0) 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(); unsigned solve_with_tableau();
bool basis_column_is_set_correctly(unsigned j) const { 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); lp_assert(this->m_basis_heading[j] >= 0);
unsigned i = static_cast<unsigned>(this->m_basis_heading[j]); unsigned i = static_cast<unsigned>(this->m_basis_heading[j]);
for (const row_cell<T> & rc : this->m_A.m_rows[i]) { for (const row_cell<T> & rc : this->m_A.m_rows[i]) {
unsigned k = rc.m_j; if (rc.dead()) continue;
unsigned k = rc.var();
if (k == j) if (k == j)
continue; continue;
this->m_d[k] += delta * rc.get_val(); this->m_d[k] += delta * rc.get_val();

2
src/util/lp/lp_primal_core_solver_def.h
View file

@ -975,7 +975,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::delete_fa
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_norms() { template <typename T, typename X> void lp_primal_core_solver<T, X>::init_column_norms() {
lp_assert(numeric_traits<T>::precise() == false); lp_assert(numeric_traits<T>::precise() == false);
for (unsigned j = 0; j < this->m_n(); j++) { for (unsigned j = 0; j < this->m_n(); j++) {
this->m_column_norms[j] = T(static_cast<int>(this->m_A.m_columns[j].size() + 1)) this->m_column_norms[j] = T(static_cast<int>(this->m_A.m_columns[j].live_size() + 1))
+ T(static_cast<int>(this->m_settings.random_next() % 10000)) / T(100000); + T(static_cast<int>(this->m_settings.random_next() % 10000)) / T(100000);
} }

22
src/util/lp/lp_primal_core_solver_tableau_def.h
View file

@ -266,17 +266,18 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
unsigned row_min_nz = this->m_n() + 1; unsigned row_min_nz = this->m_n() + 1;
m_leaving_candidates.clear(); m_leaving_candidates.clear();
auto & col = this->m_A.m_columns[entering]; 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++) { for (;k < col_size && unlimited; k++) {
const column_cell & c = col[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); const T & ed = this->m_A.get_val(c);
lp_assert(!numeric_traits<T>::is_zero(ed)); lp_assert(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i]; unsigned j = this->m_basis[i];
limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, t, unlimited); limit_theta_on_basis_column(j, - ed * m_sign_of_entering_delta, t, unlimited);
if (!unlimited) { if (!unlimited) {
m_leaving_candidates.push_back(j); 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) { if (unlimited) {
@ -288,14 +289,15 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
X ratio; X ratio;
for (;k < col_size; k++) { for (;k < col_size; k++) {
const column_cell & c = col[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); const T & ed = this->m_A.get_val(c);
lp_assert(!numeric_traits<T>::is_zero(ed)); lp_assert(!numeric_traits<T>::is_zero(ed));
unsigned j = this->m_basis[i]; unsigned j = this->m_basis[i];
unlimited = true; unlimited = true;
limit_theta_on_basis_column(j, -ed * m_sign_of_entering_delta, ratio, unlimited); limit_theta_on_basis_column(j, -ed * m_sign_of_entering_delta, ratio, unlimited);
if (unlimited) continue; 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) { if (ratio < t) {
t = ratio; t = ratio;
m_leaving_candidates.clear(); m_leaving_candidates.clear();
@ -304,7 +306,7 @@ template <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
} else if (ratio == t && i_nz < row_min_nz) { } else if (ratio == t && i_nz < row_min_nz) {
m_leaving_candidates.clear(); m_leaving_candidates.clear();
m_leaving_candidates.push_back(j); 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) { } else if (ratio == t && i_nz == row_min_nz) {
m_leaving_candidates.push_back(j); m_leaving_candidates.push_back(j);
} }
@ -348,6 +350,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
template <typename T, typename X> bool lp_primal_core_solver<T, X>:: template <typename T, typename X> bool lp_primal_core_solver<T, X>::
update_basis_and_x_tableau(int entering, int leaving, X const & tt) { update_basis_and_x_tableau(int entering, int leaving, X const & tt) {
lp_assert(this->use_tableau()); lp_assert(this->use_tableau());
lp_assert(entering != leaving);
update_x_tableau(entering, tt); update_x_tableau(entering, tt);
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
this->change_basis(entering, 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); this->add_delta_to_x_and_call_tracker(entering, delta);
if (!this->m_using_infeas_costs) { if (!this->m_using_infeas_costs) {
for (const auto & c : this->m_A.m_columns[entering]) { 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)); 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 } 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<T>()); lp_assert(this->m_costs[entering] == zero_of_type<T>());
// m_d[entering] can change because of the cost change for basic columns. // m_d[entering] can change because of the cost change for basic columns.
for (const auto & c : this->m_A.m_columns[entering]) { 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]; unsigned j = this->m_basis[i];
this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c)); this->add_delta_to_x_and_call_tracker(j, -delta * this->m_A.get_val(c));
update_inf_cost_for_column_tableau(j); update_inf_cost_for_column_tableau(j);
@ -407,7 +411,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_reduced
else { else {
T& d = this->m_d[j] = this->m_costs[j]; T& d = this->m_d[j] = this->m_costs[j];
for (auto & cc : this->m_A.m_columns[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);
} }
} }
} }

4
src/util/lp/lu.h
View file

@ -330,11 +330,11 @@ public:
for (unsigned i = m_prev; i < m; i++) { for (unsigned i = m_prev; i < m; i++) {
for (const row_cell<T> & c : m_A.m_rows[i]) { for (const row_cell<T> & c : m_A.m_rows[i]) {
int h = heading[c.m_j]; int h = heading[c.var()];
if (h < 0) { if (h < 0) {
continue; continue;
} }
columns_to_replace.insert(c.m_j); columns_to_replace.insert(c.var());
} }
} }
return columns_to_replace; return columns_to_replace;

2
src/util/lp/lu_def.h
View file

@ -386,7 +386,7 @@ void lu< M>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector<i
auto & row = m_A.m_rows[k]; auto & row = m_A.m_rows[k];
const T & y_k = y.m_data[k]; const T & y_k = y.m_data[k];
for (auto & c : row) { for (auto & c : row) {
unsigned j = c.m_j; unsigned j = c.var();
int hj = heading[j]; int hj = heading[j];
if (hj < 0) continue; if (hj < 0) continue;
if (m_ii.m_data[hj] == 0) if (m_ii.m_data[hj] == 0)

5
src/util/lp/random_updater_def.h
View file

@ -83,8 +83,9 @@ void random_updater::add_column_to_sets(unsigned j) {
add_value(m_lar_solver.get_core_solver().m_r_x[j]); add_value(m_lar_solver.get_core_solver().m_r_x[j]);
} else { } else {
unsigned row = m_lar_solver.get_core_solver().m_r_heading[j]; 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]) { for (auto & row_c : m_lar_solver.get_core_solver().m_r_A.m_rows[row].m_cells) {
unsigned cj = row_c.m_j; if (row_c.dead()) continue;
unsigned cj = row_c.var();
if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) { if (m_lar_solver.get_core_solver().m_r_heading[cj] < 0) {
m_var_set.insert(cj); m_var_set.insert(cj);
add_value(m_lar_solver.get_core_solver().m_r_x[cj]); add_value(m_lar_solver.get_core_solver().m_r_x[cj]);

9
src/util/lp/static_matrix.cpp
View file

@ -30,6 +30,9 @@ Revision History:
#include "util/lp/lar_solver.h" #include "util/lp/lar_solver.h"
namespace lp { namespace lp {
template void static_matrix<double, double>::add_columns_at_the_end(unsigned int); template void static_matrix<double, double>::add_columns_at_the_end(unsigned int);
template void static_matrix<double, double>::add_new_element(unsigned i, unsigned j, const double & v);
template void static_matrix<mpq, mpq>::add_new_element(unsigned i, unsigned j, const mpq & v);
template void static_matrix<mpq, impq>::add_new_element(unsigned i, unsigned j, const mpq & v);
template void static_matrix<double, double>::clear(); template void static_matrix<double, double>::clear();
#ifdef Z3DEBUG #ifdef Z3DEBUG
template bool static_matrix<double, double>::is_correct() const; template bool static_matrix<double, double>::is_correct() const;
@ -47,7 +50,6 @@ template double static_matrix<double, double>::get_min_abs_in_row(unsigned int)
template void static_matrix<double, double>::init_empty_matrix(unsigned int, unsigned int); template void static_matrix<double, double>::init_empty_matrix(unsigned int, unsigned int);
template void static_matrix<double, double>::init_row_columns(unsigned int, unsigned int); template void static_matrix<double, double>::init_row_columns(unsigned int, unsigned int);
template static_matrix<double, double>::ref & static_matrix<double, double>::ref::operator=(double const&); template static_matrix<double, double>::ref & static_matrix<double, double>::ref::operator=(double const&);
template void static_matrix<double, double>::set(unsigned int, unsigned int, double const&);
template static_matrix<double, double>::static_matrix(unsigned int, unsigned int); template static_matrix<double, double>::static_matrix(unsigned int, unsigned int);
template void static_matrix<mpq, mpq>::add_column_to_vector(mpq const&, unsigned int, mpq*) const; template void static_matrix<mpq, mpq>::add_column_to_vector(mpq const&, unsigned int, mpq*) const;
template void static_matrix<mpq, mpq>::add_columns_at_the_end(unsigned int); template void static_matrix<mpq, mpq>::add_columns_at_the_end(unsigned int);
@ -63,7 +65,6 @@ template mpq static_matrix<mpq, mpq>::get_min_abs_in_column(unsigned int) const;
template mpq static_matrix<mpq, mpq>::get_min_abs_in_row(unsigned int) const; template mpq static_matrix<mpq, mpq>::get_min_abs_in_row(unsigned int) const;
template void static_matrix<mpq, mpq>::init_row_columns(unsigned int, unsigned int); template void static_matrix<mpq, mpq>::init_row_columns(unsigned int, unsigned int);
template static_matrix<mpq, mpq>::ref& static_matrix<mpq, mpq>::ref::operator=(mpq const&); template static_matrix<mpq, mpq>::ref& static_matrix<mpq, mpq>::ref::operator=(mpq const&);
template void static_matrix<mpq, mpq>::set(unsigned int, unsigned int, mpq const&);
template static_matrix<mpq, mpq>::static_matrix(unsigned int, unsigned int); template static_matrix<mpq, mpq>::static_matrix(unsigned int, unsigned int);
#ifdef Z3DEBUG #ifdef Z3DEBUG
@ -72,13 +73,11 @@ template bool static_matrix<mpq, numeric_pair<mpq> >::is_correct() const;
template void static_matrix<mpq, numeric_pair<mpq> >::copy_column_to_indexed_vector(unsigned int, indexed_vector<mpq>&) const; template void static_matrix<mpq, numeric_pair<mpq> >::copy_column_to_indexed_vector(unsigned int, indexed_vector<mpq>&) const;
template mpq static_matrix<mpq, numeric_pair<mpq> >::get_elem(unsigned int, unsigned int) const; template mpq static_matrix<mpq, numeric_pair<mpq> >::get_elem(unsigned int, unsigned int) const;
template void static_matrix<mpq, numeric_pair<mpq> >::init_empty_matrix(unsigned int, unsigned int); template void static_matrix<mpq, numeric_pair<mpq> >::init_empty_matrix(unsigned int, unsigned int);
template void static_matrix<mpq, numeric_pair<mpq> >::set(unsigned int, unsigned int, mpq const&);
template bool lp::static_matrix<double, double>::pivot_row_to_row_given_cell(unsigned int, column_cell &, unsigned int); template bool lp::static_matrix<double, double>::pivot_row_to_row_given_cell(unsigned int, column_cell &, unsigned int);
template bool lp::static_matrix<lp::mpq, lp::mpq>::pivot_row_to_row_given_cell(unsigned int, column_cell& , unsigned int); template bool lp::static_matrix<lp::mpq, lp::mpq>::pivot_row_to_row_given_cell(unsigned int, column_cell& , unsigned int);
template bool lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >::pivot_row_to_row_given_cell(unsigned int, column_cell&, unsigned int); template bool lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >::pivot_row_to_row_given_cell(unsigned int, column_cell&, unsigned int);
template void lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_element(vector<lp::row_cell<lp::mpq>, true, unsigned int>&, lp::row_cell<lp::mpq>&); template void lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_element(lp::row_cell<lp::mpq>&);
} }

549
src/util/lp/static_matrix.h
View file

@ -1,22 +1,22 @@
/*++ /*++
Copyright (c) 2017 Microsoft Corporation Copyright (c) 2017 Microsoft Corporation
Module Name: Module Name:
<name> <name>
Abstract: Abstract:
<abstract> <abstract>
Author: Author:
Lev Nachmanson (levnach) Lev Nachmanson (levnach)
Revision History: Revision History:
--*/ --*/
#pragma once #pragma once
#include "util/vector.h" #include "util/vector.h"
@ -29,30 +29,387 @@ Revision History:
#include <stack> #include <stack>
namespace lp { namespace lp {
struct column_cell { class column_cell {
unsigned m_i; // points to the row unsigned m_i; // points to the row
unsigned m_offset; // the offset of the element in the matrix row unsigned m_offset; // the offset of the element in the matrix row, or the next dead cell in the column_strip
column_cell(unsigned i, unsigned offset) : m_i(i), m_offset(offset) { 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<unsigned>(-1); }
void set_dead() { m_i = -1;}
}; };
template <typename T> template <typename T>
struct row_cell { class row_cell {
unsigned m_j; // points to the column unsigned m_j; // points to the column
unsigned m_offset; // offset in column unsigned m_offset; // offset in column, or the next dead cell in the row_strip
T m_value; T m_value;
public:
row_cell(unsigned j, unsigned offset, T const & val) : m_j(j), m_offset(offset), m_value(val) { 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;} row_cell(unsigned j, T const & val) : m_j(j), m_value(val) {
T & get_val() { return m_value;} }
void set_val(const T& v) { m_value = v; } const T & get_val() const {
unsigned var() const { return m_j;} lp_assert(alive());
const T & coeff() const { return m_value;} 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<unsigned>(-1); }
void set_dead() { m_j = static_cast<unsigned>(-1); }
}; };
template <typename T> template <typename T>
using row_strip = vector<row_cell<T>>; class row_strip {
unsigned m_live_size;
int m_first_dead;
public:
row_strip() : m_live_size(0), m_first_dead(-1) {}
row_cell<T>* begin() { return m_cells.begin();}
const row_cell<T>* begin() const { return m_cells.begin();}
row_cell<T>* end() { return m_cells.end();}
const row_cell<T>* end() const { return m_cells.end();}
unsigned live_size() const { return m_live_size; }
vector<row_cell<T>> m_cells;
unsigned cells_size() const { return m_cells.size(); }
const row_cell<T> & operator[](unsigned i) const { return m_cells[i]; }
row_cell<T> & 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<int>(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<unsigned> d0;
std::set<unsigned> 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<unsigned>(-1) ; i = m_cells[i].next_dead_index())
d1.insert(i);
return d0 == d1;
}
row_cell<T> & 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<T>(j, val));
return m_cells.back();
}
void shrink_to_live() {
m_cells.shrink(live_size());
m_first_dead = -1;
}
};
class column_strip {
vector<column_cell> 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<int>(last))
m_first_dead = m_cells.back().next_dead_index();
else {
skip_dead_cell(last);
}
m_cells.pop_back();
}
std::set<unsigned> d0;
std::set<unsigned> 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<unsigned>(-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 <typename A, typename B>
void compress_cells(A& cells, vector<B>& 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!!! // each assignment for this matrix should be issued only once!!!
template <typename T, typename X> template <typename T, typename X>
@ -66,13 +423,13 @@ class static_matrix
unsigned m_n; unsigned m_n;
dim(unsigned m, unsigned n) :m_m(m), m_n(n) {} dim(unsigned m, unsigned n) :m_m(m), m_n(n) {}
}; };
std::stack<dim> m_stack; // fields
std::stack<dim> m_stack;
public: public:
typedef vector<column_cell> column_strip; vector<int> m_vector_of_row_offsets;
vector<int> m_vector_of_row_offsets; indexed_vector<T> m_work_vector;
indexed_vector<T> m_work_vector; vector<row_strip<T>> m_rows;
vector<row_strip<T>> m_rows; vector<column_strip> m_columns;
vector<column_strip> m_columns;
// starting inner classes // starting inner classes
class ref { class ref {
static_matrix & m_matrix; static_matrix & m_matrix;
@ -80,9 +437,9 @@ public:
unsigned m_col; unsigned m_col;
public: public:
ref(static_matrix & m, unsigned row, unsigned col):m_matrix(m), m_row(row), m_col(col) {} 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); } operator T () const { return m_matrix.get_elem(m_row, m_col); }
}; };
@ -98,7 +455,7 @@ public:
public: public:
const T & get_val(const column_cell & c) const { 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<T> &rc) { column_cell & get_column_cell(const row_cell<T> &rc) {
@ -107,7 +464,7 @@ public:
void init_row_columns(unsigned m, unsigned n); void init_row_columns(unsigned m, unsigned n);
// constructor with no parameters // constructor with no parameters
static_matrix() {} static_matrix() {}
// constructor // constructor
@ -142,12 +499,12 @@ public:
void remove_last_column(unsigned j); void remove_last_column(unsigned j);
void remove_element(vector<row_cell<T>> & row, row_cell<T> & elem_to_remove); void remove_element(row_cell<T> & elem_to_remove);
void multiply_column(unsigned column, T const & alpha) { void multiply_column(unsigned column, T const & alpha) {
for (auto & t : m_columns[column]) { for (auto & t : m_columns[column]) {
auto & r = m_rows[t.m_i][t.m_offset]; auto & r = m_rows[t.var()].m_cells[t.offset()];
r.m_value *= alpha; r.coeff() *= alpha;
} }
} }
@ -165,8 +522,6 @@ public:
return column_cell(column, offset); 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); } ref operator()(unsigned row, unsigned col) { return ref(*this, row, col); }
std::set<std::pair<unsigned, unsigned>> get_domain(); std::set<std::pair<unsigned, unsigned>> get_domain();
@ -175,8 +530,8 @@ public:
T get_max_abs_in_row(unsigned row) const; T get_max_abs_in_row(unsigned row) const;
void add_column_to_vector (const T & a, unsigned j, T * v) const { void add_column_to_vector (const T & a, unsigned j, T * v) const {
for (const auto & it : m_columns[j]) { for (const auto & c : m_columns[j]) {
v[it.m_i] += a * get_val(it); v[c.var()] += a * get_val(c);
} }
} }
@ -189,9 +544,6 @@ public:
void check_consistency(); void check_consistency();
#endif #endif
void cross_out_row(unsigned k);
// //
void fix_row_indices_in_each_column_for_crossed_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; 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 number_of_non_zeroes() const {
unsigned ret = 0; unsigned ret = 0;
@ -237,12 +589,12 @@ public:
m_stack.push(d); m_stack.push(d);
} }
void pop_row_columns(const vector<row_cell<T>> & row) { void pop_row_columns(const row_strip<T> & row) {
for (auto & c : row) { for (auto & c : row) {
unsigned j = c.m_j; if (c.dead()) {
auto & col = m_columns[j]; continue;
lp_assert(col[col.size() - 1].m_i == m_rows.size() -1 ); // todo : start here!!!! }
col.pop_back(); m_columns[c.var()].delete_at(c.offset());
} }
} }
@ -272,8 +624,9 @@ public:
} }
void multiply_row(unsigned row, T const & alpha) { void multiply_row(unsigned row, T const & alpha) {
for (auto & t : m_rows[row]) { for (auto & t : m_rows[row].m_cells) {
t.m_value *= alpha; if (t.alive())
t.coeff() *= alpha;
} }
} }
@ -286,8 +639,8 @@ public:
T dot_product_with_column(const vector<T> & y, unsigned j) const { T dot_product_with_column(const vector<T> & y, unsigned j) const {
lp_assert(j < column_count()); lp_assert(j < column_count());
T ret = numeric_traits<T>::zero(); T ret = numeric_traits<T>::zero();
for (auto & it : m_columns[j]) { for (auto & c : m_columns[j]) {
ret += y[it.m_i] * get_val(it); // get_value_of_column_cell(it); ret += y[c.var()] * get_val(c);
} }
return ret; return ret;
} }
@ -301,18 +654,20 @@ public:
m_rows[i] = m_rows[ii]; m_rows[i] = m_rows[ii];
m_rows[ii] = t; m_rows[ii] = t;
// now fix the columns // now fix the columns
for (auto & rc : m_rows[i]) { for (const auto & rc : m_rows[i]) {
column_cell & cc = m_columns[rc.m_j][rc.m_offset]; if (rc.dead()) continue;
lp_assert(cc.m_i == ii); column_cell & cc = m_columns[rc.var()][rc.offset()];
cc.m_i = i; lp_assert(cc.var() == ii);
cc.var() = i;
} }
for (auto & rc : m_rows[ii]) { for (const auto & rc : m_rows[ii]) {
column_cell & cc = m_columns[rc.m_j][rc.m_offset]; if (rc.dead()) continue;
lp_assert(cc.m_i == i); column_cell & cc = m_columns[rc.var()][rc.offset()];
cc.m_i = ii; lp_assert(cc.var() == i);
cc.var() = ii;
} }
} }
void fill_last_row_with_pivoting_loop_block(unsigned j, const vector<int> & basis_heading) { void fill_last_row_with_pivoting_loop_block(unsigned j, const vector<int> & basis_heading) {
int row_index = basis_heading[j]; int row_index = basis_heading[j];
if (row_index < 0) if (row_index < 0)
@ -322,17 +677,18 @@ public:
return; return;
for (const auto & c : m_rows[row_index]) { for (const auto & c : m_rows[row_index]) {
if (c.m_j == j) { if (c.dead()) continue;
if (c.var() == j) {
continue; 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); bool was_zero = is_zero(wv);
wv -= alpha * c.m_value; wv -= alpha * c.coeff();
if (was_zero) if (was_zero)
m_work_vector.m_index.push_back(c.m_j); m_work_vector.m_index.push_back(c.var());
else { else {
if (is_zero(wv)) { 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); lp_assert(row_count() > 0);
m_work_vector.resize(column_count()); m_work_vector.resize(column_count());
T a; T a;
// we use the form -it + 1 = 0 // we use the form -it + 1 = 0
m_work_vector.set_value(one_of_type<T>(), bj); m_work_vector.set_value(one_of_type<T>(), bj);
for (auto p : row) { for (auto p : row) {
m_work_vector.set_value(-p.coeff(), p.var()); m_work_vector.set_value(-p.coeff(), p.var());
@ -366,10 +722,10 @@ public:
unsigned last_row = row_count() - 1; unsigned last_row = row_count() - 1;
for (unsigned j : m_work_vector.m_index) { 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); lp_assert(column_count() > 0);
set(last_row, column_count() - 1, one_of_type<T>()); add_new_element(last_row, column_count() - 1, one_of_type<T>());
} }
void copy_column_to_vector (unsigned j, vector<T> & v) const { void copy_column_to_vector (unsigned j, vector<T> & v) const {
@ -377,7 +733,7 @@ public:
for (auto & it : m_columns[j]) { for (auto & it : m_columns[j]) {
const T& val = get_val(it); const T& val = get_val(it);
if (!is_zero(val)) if (!is_zero(val))
v[it.m_i] = val; v[it.var()] = val;
} }
} }
@ -386,7 +742,8 @@ public:
L ret = zero_of_type<L>(); L ret = zero_of_type<L>();
lp_assert(row < m_rows.size()); lp_assert(row < m_rows.size());
for (auto & it : m_rows[row]) { 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; return ret;
} }
@ -397,9 +754,9 @@ public:
// constructor // constructor
column_cell_plus(const column_cell & c, const static_matrix& A) : column_cell_plus(const column_cell & c, const static_matrix& A) :
m_c(c), m_A(A) {} m_c(c), m_A(A) {}
unsigned var() const { return m_c.m_i; } unsigned var() const { return m_c.var(); }
const T & coeff() const { return m_A.m_rows[var()][m_c.m_offset].get_val(); } 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 { struct column_container {
@ -411,7 +768,7 @@ public:
// fields // fields
const column_cell *m_c; const column_cell *m_c;
const static_matrix& m_A; const static_matrix& m_A;
const column_cell *m_end;
//typedefs //typedefs
@ -425,13 +782,19 @@ public:
reference operator*() const { reference operator*() const {
return column_cell_plus(*m_c, m_A); return column_cell_plus(*m_c, m_A);
} }
self_type operator++() { self_type i = *this; m_c++; return i; } self_type operator++() { self_type i = *this;
self_type operator++(int) { m_c++; return *this; } m_c++;
return i;
}
self_type operator++(int) {
m_c++;
return *this;
}
const_iterator(const column_cell* it, const static_matrix& A) : const_iterator(const column_cell* it, const static_matrix& A) :
m_c(it), m_c(it),
m_A(A) m_A(A){}
{}
bool operator==(const self_type &other) const { bool operator==(const self_type &other) const {
return m_c == other.m_c; return m_c == other.m_c;
} }
@ -449,8 +812,32 @@ public:
column_container column(unsigned j) const { column_container column(unsigned j) const {
return column_container(j, *this); 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);} ref_row operator[](unsigned i) const { return ref_row(*this, i);}
typedef T coefftype; typedef T coefftype;
typedef X argtype; typedef X argtype;

204
src/util/lp/static_matrix_def.h
View file

@ -36,44 +36,54 @@ void static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(const row_strip<T> & rvals) { template <typename T, typename X> void static_matrix<T, X>::scan_row_ii_to_offset_vector(const row_strip<T> & rvals) {
for (unsigned j = 0; j < rvals.size(); j++) for (unsigned j = 0; j < rvals.cells_size(); j++) {
m_vector_of_row_offsets[rvals[j].m_j] = j; if (rvals[j].dead()) continue;
m_vector_of_row_offsets[rvals[j].var()] = j;
}
} }
template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_given_cell(unsigned i, column_cell & c, unsigned pivot_col) { template <typename T, typename X> bool static_matrix<T, X>::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); lp_assert(i < row_count() && ii < column_count() && i != ii);
T alpha = -get_val(c); T alpha = -get_val(c);
lp_assert(!is_zero(alpha)); lp_assert(!is_zero(alpha));
compress_row_if_needed(ii);
auto & rowii = m_rows[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); 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 // run over the pivot row and update row ii
for (const auto & iv : m_rows[i]) { for (const auto & iv : m_rows[i].m_cells) {
unsigned j = iv.m_j; if (iv.dead()) continue;
unsigned j = iv.var();
if (j == pivot_col) continue; if (j == pivot_col) continue;
T alv = alpha * iv.m_value; lp_assert(!is_zero(iv.coeff()));
lp_assert(!is_zero(iv.m_value)); T alv = alpha * iv.coeff();
int j_offs = m_vector_of_row_offsets[j]; int j_offs = m_vector_of_row_offsets[j];
if (j_offs == -1) { // it is a new element if (j_offs == -1) { // it is a new element
add_new_element(ii, j, alv); add_new_element(ii, j, alv);
} }
else { else {
rowii[j_offs].m_value += alv; rowii[j_offs].coeff() += alv;
} }
} }
// clean the work vector // clean the work vector
for (unsigned k = 0; k < prev_size_ii; k++) { 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 // remove zeroes
for (unsigned k = rowii.size(); k-- > 0; ) { for (unsigned k = rowii.cells_size(); k-- > 0; ) {
if (is_zero(rowii[k].m_value)) auto & c = rowii[k];
remove_element(rowii, rowii[k]); if (c.dead())
continue;
if (is_zero(c.coeff()))
remove_element(c);
} }
lp_assert(is_correct());
return !rowii.empty(); return !rowii.empty();
} }
@ -86,7 +96,7 @@ static_matrix<T, X>::static_matrix(static_matrix const &A, unsigned * /* basis *
init_row_columns(m, m); init_row_columns(m, m);
while (m--) { while (m--) {
for (auto & col : A.m_columns[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 <typename T, typename X> void static_matrix<T, X>::init_empty_matrix
init_row_columns(m, n); init_row_columns(m, n);
} }
template <typename T, typename X> unsigned static_matrix<T, X>::lowest_row_in_column(unsigned col) { template <typename T, typename X> unsigned static_matrix<T, X>::lowest_row_in_column(unsigned col) {
lp_assert(col < column_count());
column_strip & colstrip = m_columns[col]; column_strip & colstrip = m_columns[col];
lp_assert(colstrip.size() > 0); lp_assert(colstrip.live_size() > 0);
unsigned ret = 0; unsigned ret = 0;
for (auto & t : colstrip) { for (const auto & t : colstrip) {
if (t.m_i > ret) { if (t.dead()) continue;
ret = t.m_i; if (t.var() > ret) {
ret = t.var();
} }
} }
return ret; return ret;
@ -136,10 +146,10 @@ template <typename T, typename X> void static_matrix<T, X>::forget_last_colum
template <typename T, typename X> void static_matrix<T, X>::remove_last_column(unsigned j) { template <typename T, typename X> void static_matrix<T, X>::remove_last_column(unsigned j) {
column_strip & col = m_columns.back(); column_strip & col = m_columns.back();
for (auto & it : col) { for (auto & it : col) {
auto & row = m_rows[it.m_i]; auto & row = m_rows[it.var()];
unsigned offset = row.size() - 1; unsigned offset = row.size() - 1;
for (auto row_it = row.rbegin(); row_it != row.rend(); row_it ++) { 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); row.erase(row.begin() + offset);
break; break;
} }
@ -150,24 +160,12 @@ template <typename T, typename X> void static_matrix<T, X>::remove_last_column(u
m_vector_of_row_offsets.pop_back(); m_vector_of_row_offsets.pop_back();
} }
template <typename T, typename X> void static_matrix<T, X>::set(unsigned row, unsigned col, T const & val) {
if (numeric_traits<T>::is_zero(val)) return;
lp_assert(row < row_count() && col < column_count());
auto & r = m_rows[row];
unsigned offs_in_cols = static_cast<unsigned>(m_columns[col].size());
m_columns[col].push_back(make_column_cell(row, static_cast<unsigned>(r.size())));
r.push_back(make_row_cell(col, offs_in_cols, val));
}
template <typename T, typename X> template <typename T, typename X>
std::set<std::pair<unsigned, unsigned>> static_matrix<T, X>::get_domain() { std::set<std::pair<unsigned, unsigned>> static_matrix<T, X>::get_domain() {
std::set<std::pair<unsigned, unsigned>> ret; std::set<std::pair<unsigned, unsigned>> ret;
for (unsigned i = 0; i < m_rows.size(); i++) { for (unsigned i = 0; i < m_rows.size(); i++) {
for (auto &it : m_rows[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; return ret;
@ -179,7 +177,7 @@ template <typename T, typename X> void static_matrix<T, X>::copy_column_to_in
for (auto & it : m_columns[j]) { for (auto & it : m_columns[j]) {
const T& val = get_val(it); const T& val = get_val(it);
if (!is_zero(val)) if (!is_zero(val))
v.set_value(val, it.m_i); v.set_value(val, it.var());
} }
} }
@ -243,7 +241,7 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency
std::unordered_map<std::pair<unsigned, unsigned>, T> by_rows; std::unordered_map<std::pair<unsigned, unsigned>, T> by_rows;
for (int i = 0; i < m_rows.size(); i++){ for (int i = 0; i < m_rows.size(); i++){
for (auto & t : m_rows[i]) { for (auto & t : m_rows[i]) {
std::pair<unsigned, unsigned> p(i, t.m_j); std::pair<unsigned, unsigned> p(i, t.var());
lp_assert(by_rows.find(p) == by_rows.end()); lp_assert(by_rows.find(p) == by_rows.end());
by_rows[p] = t.get_val(); by_rows[p] = t.get_val();
} }
@ -251,7 +249,7 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency
std::unordered_map<std::pair<unsigned, unsigned>, T> by_cols; std::unordered_map<std::pair<unsigned, unsigned>, T> by_cols;
for (int i = 0; i < m_columns.size(); i++){ for (int i = 0; i < m_columns.size(); i++){
for (auto & t : m_columns[i]) { for (auto & t : m_columns[i]) {
std::pair<unsigned, unsigned> p(t.m_i, i); std::pair<unsigned, unsigned> p(t.var(), i);
lp_assert(by_cols.find(p) == by_cols.end()); lp_assert(by_cols.find(p) == by_cols.end());
by_cols[p] = get_val(t); by_cols[p] = get_val(t);
} }
@ -266,51 +264,10 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency
} }
#endif #endif
template <typename T, typename X> void static_matrix<T, X>::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 <typename T, typename X> void static_matrix<T, X>::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 <typename T, typename X> void static_matrix<T, X>::cross_out_row_from_columns(unsigned k, row_strip<T> & row) {
for (auto & t : row) {
cross_out_row_from_column(t.m_j, k);
}
}
template <typename T, typename X> void static_matrix<T, X>::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 <typename T, typename X> T static_matrix<T, X>::get_elem(unsigned i, unsigned j) const { // should not be used in efficient code !!!! template <typename T, typename X> T static_matrix<T, X>::get_elem(unsigned i, unsigned j) const { // should not be used in efficient code !!!!
for (auto & t : m_rows[i]) { for (auto & t : m_rows[i]) {
if (t.m_j == j) { if (t.dead()) continue;
if (t.var() == j) {
return t.get_val(); return t.get_val();
} }
} }
@ -342,16 +299,22 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
auto &r = m_rows[i]; auto &r = m_rows[i];
std::unordered_set<unsigned> s; std::unordered_set<unsigned> s;
for (auto & rc : r) { 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; return false;
} }
s.insert(rc.m_j); s.insert(rc.var());
if (rc.m_j >= m_columns.size()) if (rc.var() >= m_columns.size())
return false; 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; 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; return false;
if (& m_rows[cc.var()][cc.offset()] != & rc) {
return false;
}
if (is_zero(rc.get_val())) { if (is_zero(rc.get_val())) {
return false; return false;
} }
@ -363,51 +326,56 @@ template <typename T, typename X> bool static_matrix<T, X>::is_correct() const {
auto & c = m_columns[j]; auto & c = m_columns[j];
std::unordered_set<unsigned> s; std::unordered_set<unsigned> s;
for (auto & cc : c) { 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; return false;
} }
s.insert(cc.m_i); s.insert(cc.var());
if (cc.m_i >= m_rows.size()) if (cc.var() >= m_rows.size())
return false; 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; 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; 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; return true;
} }
template <typename T, typename X> template <typename T, typename X>
void static_matrix<T, X>::remove_element(vector<row_cell<T>> & row_vals, row_cell<T> & row_el_iv) { void static_matrix<T, X>::remove_element(row_cell<T> & rc) {
unsigned column_offset = row_el_iv.m_offset; lp_assert(rc.alive());
auto & column_vals = m_columns[row_el_iv.m_j]; unsigned j = rc.var();
column_cell& cs = m_columns[row_el_iv.m_j][column_offset]; unsigned j_offset = rc.offset();
unsigned row_offset = cs.m_offset; auto & col = m_columns[j];
if (column_offset != column_vals.size() - 1) { column_cell & c = col[j_offset];
auto & cc = column_vals[column_offset] = column_vals.back(); // copy from the tail unsigned i = c.var();
m_rows[cc.m_i][cc.m_offset].m_offset = column_offset; unsigned i_offset = c.offset();
} col.delete_at(j_offset);
m_rows[i].delete_at(i_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();
} }
template <typename T, typename X> template <typename T, typename X>
void static_matrix<T, X>::add_new_element(unsigned row, unsigned col, const T& val) { void static_matrix<T, X>::add_new_element(unsigned i, unsigned j, const T& val) {
auto & row_vals = m_rows[row]; auto & row = m_rows[i];
auto & col_vals = m_columns[col]; auto & col = m_columns[j];
unsigned row_el_offs = static_cast<unsigned>(row_vals.size()); unsigned offset_in_row, offset_in_col;
unsigned col_el_offs = static_cast<unsigned>(col_vals.size()); row_cell<T>& rc = row.add_cell(j, val, offset_in_row);
row_vals.push_back(row_cell<T>(col, col_el_offs, val)); column_cell& cc = col.add_cell(i, offset_in_col);
col_vals.push_back(column_cell(row, row_el_offs)); rc.offset() = offset_in_col;
cc.offset() = offset_in_row;
} }
} }