mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-12 12:08:18 +00:00
Code Activity

dio with static_matrix initial setup

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2024-10-10 19:23:37 -07:00 committed by Lev Nachmanson
parent 9e8b17b5f8
commit 42bdc893a9
4 changed files with 363 additions and 243 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

513
src/math/lp/dioph_eq.cpp
View file

@ -56,6 +56,7 @@ namespace lp {
} }
}; };
std::ostream& print_S(std::ostream & out) { std::ostream& print_S(std::ostream & out) {
out << "S:\n"; out << "S:\n";
for (unsigned i : m_s) { for (unsigned i : m_s) {
@ -66,7 +67,8 @@ namespace lp {
} }
std::ostream& print_lar_term_L(const lar_term & t, std::ostream & out) const { std::ostream& print_lar_term_L(const lar_term & t, std::ostream & out) const {
return print_linear_combination_customized(t.coeffs_as_vector(), [](int j)->std::string {return "y"+std::to_string(j);}, out ); return print_linear_combination_customized(t.coeffs_as_vector(),
[](int j)->std::string {return "y"+std::to_string(j);}, out );
} }
std::ostream& print_term_o(term_o const& term, std::ostream& out) const { std::ostream& print_term_o(term_o const& term, std::ostream& out) const {
@ -110,26 +112,29 @@ namespace lp {
An annotated state is a triple ⟨E′, λ, σ⟩, where E′ is a set of pairs ⟨e, ℓ⟩ in which An annotated state is a triple ⟨E′, λ, σ⟩, where E′ is a set of pairs ⟨e, ℓ⟩ in which
e is an equation and ℓ is a linear combination of variables from L e is an equation and ℓ is a linear combination of variables from L
*/ */
struct eprime_pair { //
term_o m_e; struct eprime_entry {
unsigned m_row_index; // the index of the row in the constraint matrix that m_e corresponds to
term_o m_e; // it will be used for debugging only
// we keep the dependency of the equation in m_l // we keep the dependency of the equation in m_l
// a more expensive alternative is to keep the history term of m_e : originally m_l is i, the index of row m_e was constructed from // a more expensive alternative is to keep the history term of m_e : originally m_l is i, the index of row m_e was constructed from
u_dependency *m_l; u_dependency *m_l;
}; };
vector<eprime_pair> m_eprime; enum class row_status {
F,
/* let σ be a partial mapping from variables in L united with X to linear combinations S,
of variables in X and of integer constants showing the substitutions NO_S_NO_F
*/ };
u_map<term_o> m_sigma; vector<eprime_entry> m_eprime;
// the terms are stored in m_A and m_c
public: static_matrix<mpq, numeric_pair<mpq>> m_e_matrix; // the rows of the matrix are the terms, without the constant part
vector<row_status> m_row_status;
vector<mpq> m_c; // to keep the constants of the terms
int_solver& lia; int_solver& lia;
lar_solver& lra; lar_solver& lra;
explanation m_infeas_explanation; explanation m_infeas_explanation;
indexed_vector<mpq> m_indexed_work_vector;
// we start assigning with UINT_MAX and go down, print it as l(UINT_MAX - m_last_fresh_x_var)
unsigned m_last_fresh_x_var;
bool m_report_branch = false; bool m_report_branch = false;
// set F // set F
@ -140,10 +145,20 @@ namespace lp {
// gives the order of substitution // gives the order of substitution
unsigned m_conflict_index = -1; // m_eprime[m_conflict_index] gives the conflict unsigned m_conflict_index = -1; // m_eprime[m_conflict_index] gives the conflict
public:
imp(int_solver& lia, lar_solver& lra): lia(lia), lra(lra) {} imp(int_solver& lia, lar_solver& lra): lia(lia), lra(lra) {}
term_o get_term_from_e_matrix(unsigned i) {
term_o row_to_term(const row_strip<mpq>& row) const { term_o t;
for (const auto & p: m_e_matrix.m_rows[i]) {
t.add_monomial(p.coeff(), p.var());
}
t.c() = m_c[i];
return t;
}
private:
term_o create_eprime_entry_from_row(const row_strip<mpq>& row, unsigned row_index) {
const auto lcm = get_denominators_lcm(row); const auto lcm = get_denominators_lcm(row);
#if Z3DEBUG
term_o t; term_o t;
for (const auto & p: row) for (const auto & p: row)
if (lia.is_fixed(p.var())) if (lia.is_fixed(p.var()))
@ -151,11 +166,23 @@ namespace lp {
else else
t.add_monomial(lcm * p.coeff(), p.var()); t.add_monomial(lcm * p.coeff(), p.var());
t.c() *= lcm; t.c() *= lcm;
#endif
// init m_e_matrix and m_c
mpq & c = m_c[row_index];
for (const auto & p: row)
if (lia.is_fixed(p.var()))
c += p.coeff()*lia.lower_bound(p.var()).x;
else
m_e_matrix.add_new_element(row_index, p.var(), lcm * p.coeff());
c *= lcm;
SASSERT(t == get_term_from_e_matrix(row_index));
return t; return t;
} }
void init() { void init() {
m_last_fresh_x_var = lra.column_count(); m_e_matrix = static_matrix<mpq, impq>(lra.row_count(), lra.column_count());
m_row_status.resize(lra.row_count(), row_status::NO_S_NO_F);
m_c.resize(lra.row_count(), mpq(0));
m_report_branch = false; m_report_branch = false;
unsigned n_of_rows = lra.A_r().row_count(); unsigned n_of_rows = lra.A_r().row_count();
m_k2s.clear(); m_k2s.clear();
@ -181,17 +208,19 @@ namespace lp {
TRACE("dioph_eq", tout << "not all vars are int and small\n";); TRACE("dioph_eq", tout << "not all vars are int and small\n";);
continue; continue;
} }
term_o t = row_to_term(row); term_o t = create_eprime_entry_from_row(row, i);
m_row_status[i] = row_status::F;
TRACE("dioph_eq", tout << "row = "; lra.print_row(row, tout) << std::endl;); TRACE("dioph_eq", tout << "row = "; lra.print_row(row, tout) << std::endl;);
if (t.size() == 0) { if (t.size() == 0) {
SASSERT(t.c().is_zero()); SASSERT(t.c().is_zero());
continue; continue;
} }
// eprime_pair pair = {t, lar_term(i)}; // eprime_pair pair = {t, lar_term(i)};
eprime_pair pair = {t, get_dep_from_row(row)}; eprime_entry pair = {i, t, get_dep_from_row(row)};
m_f.push_back(static_cast<unsigned>(m_f.size()));
m_f.push_back(static_cast<unsigned>(i));
m_eprime.push_back(pair); m_eprime.push_back(pair);
TRACE("dioph_eq", print_eprime_entry(static_cast<unsigned>(m_f.size()) - 1, tout);); TRACE("dioph_eq", print_eprime_entry(static_cast<unsigned>(i), tout););
} }
} }
@ -207,9 +236,9 @@ namespace lp {
return dep; return dep;
} }
mpq gcd_of_coeffs(const term_o& t) { mpq gcd_of_row(unsigned row_index) {
mpq g(0); mpq g(0);
for (const auto & p : t) { for (const auto & p : m_e_matrix.m_rows[row_index]) {
if (g.is_zero()) if (g.is_zero())
g = abs(p.coeff()); g = abs(p.coeff());
else else
@ -235,14 +264,14 @@ namespace lp {
return false; return false;
} }
void prepare_lia_branch_report(const term_o & e, const mpq& g, const mpq new_c) { void prepare_lia_branch_report(const eprime_entry & e, const mpq& g, const mpq new_c) {
/* We have ep.m_e/g = 0, or sum((coff_i/g)*x_i) + new_c = 0, /* We have ep.m_e/g = 0, or sum((coff_i/g)*x_i) + new_c = 0,
or sum((coeff_i/g)*x_i) = -new_c, where new_c is not an integer or sum((coeff_i/g)*x_i) = -new_c, where new_c is not an integer
Then sum((coeff_i/g)*x_i) <= floor(-new_c) or sum((coeff_i/g)*x_i) >= ceil(-new_c) Then sum((coeff_i/g)*x_i) <= floor(-new_c) or sum((coeff_i/g)*x_i) >= ceil(-new_c)
*/ */
lar_term& t = lia.get_term(); lar_term& t = lia.get_term();
for (const auto& p: e) { for (const auto& p: m_e_matrix.m_rows[e.m_row_index]) {
t.add_monomial(p.coeff()/g, p.j()); t.add_monomial(p.coeff()/g, p.var());
} }
lia.offset() = floor(-new_c); lia.offset() = floor(-new_c);
lia.is_upper() = true; lia.is_upper() = true;
@ -251,23 +280,24 @@ namespace lp {
} }
// returns true if no conflict is found and false otherwise // returns true if no conflict is found and false otherwise
// this function devides all cooficients, and the free constant, of the ep.m_e by the gcd of all coefficients, // this function devides all cooficients, and the free constant, of the ep.m_e by the gcd of all coefficients,
// it is needed by the next steps // it is needed by the next steps
// the conflict can be used to report "cuts from proofs" // the conflict can be used to report "cuts from proofs"
bool normalize_e_by_gcd(eprime_pair& ep) { bool normalize_e_by_gcd(unsigned row_index) {
eprime_entry& ep = m_eprime[row_index];
TRACE("dioph_eq", print_eprime_entry(ep, tout) << std::endl;); TRACE("dioph_eq", print_eprime_entry(ep, tout) << std::endl;);
mpq g = gcd_of_coeffs(ep.m_e); mpq g = gcd_of_row(row_index);
if (g.is_zero() || g.is_one()) { if (g.is_zero() || g.is_one()) {
SASSERT(g.is_one() || ep.m_e.c().is_zero()); SASSERT(g.is_one() || ep.m_e.c().is_zero());
return true; return true;
} }
TRACE("dioph_eq", tout << "g:" << g << std::endl;); TRACE("dioph_eq", tout << "g:" << g << std::endl;);
mpq c_g = ep.m_e.c() / g; mpq c_g = m_c[row_index] / g;
if (c_g.is_int()) { if (c_g.is_int()) {
for (auto& p: ep.m_e.coeffs()) { for (auto& p: m_e_matrix.m_rows[row_index]) {
p.m_value /= g; p.coeff() /= g;
} }
ep.m_e.c() = c_g; m_c[row_index] = c_g;
// ep.m_l *= (1/g); // ep.m_l *= (1/g);
TRACE("dioph_eq", tout << "ep_m_e:"; print_eprime_entry(ep, tout) << std::endl;); TRACE("dioph_eq", tout << "ep_m_e:"; print_eprime_entry(ep, tout) << std::endl;);
return true; return true;
@ -275,7 +305,7 @@ namespace lp {
// c_g is not integral // c_g is not integral
if (lra.settings().stats().m_dio_conflicts % lra.settings().dio_cut_from_proof_period() == 0 && if (lra.settings().stats().m_dio_conflicts % lra.settings().dio_cut_from_proof_period() == 0 &&
!has_fresh_var(ep.m_e)) !has_fresh_var(ep.m_e))
prepare_lia_branch_report(ep.m_e, g, c_g); prepare_lia_branch_report(ep, g, c_g);
return false; return false;
} }
@ -283,7 +313,7 @@ namespace lp {
// returns true if no conflict is found and false otherwise // returns true if no conflict is found and false otherwise
bool normalize_by_gcd() { bool normalize_by_gcd() {
for (unsigned l: m_f) { for (unsigned l: m_f) {
if (!normalize_e_by_gcd(m_eprime[l])) { if (!normalize_e_by_gcd(l)) {
m_conflict_index = l; m_conflict_index = l;
return false; return false;
} }
@ -301,58 +331,59 @@ namespace lp {
// We look at term e.m_e: it is in form (+-)x_k + sum {a_i*x_i} + c = 0. // We look at term e.m_e: it is in form (+-)x_k + sum {a_i*x_i} + c = 0.
// We substitute x_k in t by (+-)coeff*(sum {a_i*x_i} + c), where coeff is the coefficient of x_k in t. // We substitute x_k in t by (+-)coeff*(sum {a_i*x_i} + c), where coeff is the coefficient of x_k in t.
void substitute_k_with_S_entry_for_tigthening(const eprime_pair& e, unsigned k, term_o& t, std::queue<unsigned> & q) { void substitute_k_with_S_entry_for_tigthening(const eprime_entry& e, unsigned k, std::queue<unsigned> & q) {
SASSERT (t.contains(k) && e.m_e.contains(k)); // SASSERT ( e.m_e.contains(k));
mpq coeff = t.get_coeff(k); // need to copy it because the pointer value can be changed during the next loop // mpq coeff = m_indexed_work_vector[k]; // need to copy it because the pointer value can be changed during the next loop
const mpq& k_coeff_in_e = e.m_e.get_coeff(k); // const mpq& k_coeff_in_e = e.m_e.get_coeff(k); // get it from m_e_matrix instead of e.m_e
bool is_one = k_coeff_in_e.is_one();
SASSERT(is_one || k_coeff_in_e.is_minus_one()); // bool is_one = k_coeff_in_e.is_one();
t.erase_var(k); // SASSERT(is_one || k_coeff_in_e.is_minus_one());
if (is_one) { // t.erase_var(k);
coeff = -coeff; // if (is_one) {
} // coeff = -coeff;
for (const auto& p: e.m_e) { // }
if (p.j() == k) continue; // for (const auto& p: e.m_e) {
t.add_monomial(coeff*p.coeff(), p.j()); // if (p.j() == k) continue;
} // t.add_monomial(coeff*p.coeff(), p.j());
t.c() += coeff*e.m_e.c(); // }
TRACE("dioph_eq", tout << "after subs_k\n"; print_term_o(t, tout) << std::endl;); // t.c() += coeff*e.m_e.c();
for (const auto& p: t) { // TRACE("dioph_eq", tout << "after subs_k\n"; print_term_o(t, tout) << std::endl;);
if (is_fresh_var(p.j())) { // for (const auto& p: t) {
continue; // if (is_fresh_var(p.j())) {
} // continue;
if (m_k2s[p.j()] != null_lpvar) // }
q.push(p.j()); // if (m_k2s[p.j()] != null_lpvar)
} // q.push(p.j());
// }
} }
const eprime_pair& k_th_entry(unsigned k) const { const eprime_entry& k_th_entry(unsigned k) const {
return m_eprime[m_k2s[k]]; return m_eprime[m_k2s[k]];
} }
const unsigned sub_index(unsigned k) const { const unsigned sub_index(unsigned k) const {
return m_k2s[k]; return m_k2s[k];
} }
// works on m_indexed_work_vector
void substitude_term_on_q_with_S_for_tightening(std::queue<unsigned> &q, term_o& t, u_dependency* &dep) { void substitude_term_on_q_with_S_for_tightening(std::queue<unsigned> &q, u_dependency* &dep) {
TRACE("dioph_eq_dep", tout << "dep:"; print_dep(tout, dep) << std::endl;); // TRACE("dioph_eq_dep", tout << "dep:"; print_dep(tout, dep) << std::endl;);
while (!q.empty()) { // while (!q.empty()) {
unsigned k = q.front(); // unsigned k = q.front();
q.pop(); // q.pop();
const eprime_pair& e = k_th_entry(k); // const eprime_entry& e = k_th_entry(k);
if (!t.contains(k)) { // if (m_indexed_work_vector[k].is_zero()) {
continue; // continue;
} // }
TRACE("dioph_eq", tout << "k:" << k << " in: "; print_eprime_entry(sub_index(k), tout) << std::endl;); // TRACE("dioph_eq", tout << "k:" << k << " in: "; print_eprime_entry(sub_index(k), tout) << std::endl;);
substitute_k_with_S_entry_for_tigthening(e, k, t, q); // substitute_k_with_S_entry_for_tigthening(e, k, q);
TRACE("dioph_eq", print_queue(q, tout);); // TRACE("dioph_eq", print_queue(q, tout););
dep = lra.mk_join(dep, e.m_l); // dep = lra.mk_join(dep, e.m_l);
TRACE("dioph_eq", tout << "substituted t:"; print_term_o(t, tout) << std::endl; // TRACE("dioph_eq", tout << "substituted t:"; print_term_o(t, tout) << std::endl;
tout << "\ndep:"; print_dep(tout, dep) << std::endl;); // tout << "\ndep:"; print_dep(tout, dep) << std::endl;);
} // }
} }
@ -368,6 +399,7 @@ namespace lp {
if (!m_infeas_explanation.empty()) { if (!m_infeas_explanation.empty()) {
return lia_move::conflict; return lia_move::conflict;
} }
} }
if (!change) if (!change)
return lia_move::undef; return lia_move::undef;
@ -389,38 +421,38 @@ namespace lp {
return out; return out;
} }
// j is the index of the column representing a term // j is the index of the column representing a term
// return true if there is a change // return true if a new tighter bound was set on j
bool tighten_bounds_for_column(unsigned j) { bool tighten_bounds_for_column(unsigned j) {
TRACE("dioph_eq", tout << "j:"; tout << j << std::endl;); // TRACE("dioph_eq", tout << "j:"; tout << j << std::endl;);
const lar_term& lar_t = lra.get_term(j); // const lar_term& lar_t = lra.get_term(j);
TRACE("dioph_eq", tout << "tighten_term_for_S: "; print_lar_term_L(lar_t, tout) << std::endl;); // TRACE("dioph_eq", tout << "tighten_term_for_S: "; print_lar_term_L(lar_t, tout) << std::endl;);
// create a copy: t is a copy of lar_t // // create a copy: t is a copy of lar_t
term_o t; // std::queue<unsigned> q;
std::queue<unsigned> q; // m_indexed_work_vector.clear();
for (const auto& p: lar_t) { // for (const auto& p: lar_t) {
if (sub_index(p.j()) != null_lpvar) // if (sub_index(p.j()) != null_lpvar)
q.push(p.j()); // q.push(p.j());
t.add_monomial(p.coeff(), p.j()); // m_indexed_work_vector.set_value(p.coeff(), p.j());
} // }
u_dependency * dep = nullptr; // u_dependency * dep = nullptr;
TRACE("dioph_eq", tout << "t:"; print_term_o(t, tout) << std::endl; // TRACE("dioph_eq", tout << "t:"; print_term_o(t, tout) << std::endl;
/*tout << "dep:"; print_dep(tout, dep) << std::endl;*/); // /*tout << "dep:"; print_dep(tout, dep) << std::endl;*/);
substitude_term_on_q_with_S_for_tightening(q, t, dep); // substitude_term_on_q_with_S_for_tightening(q, dep);
TRACE("dioph_eq", tout << "after process_q_with_S\n t:"; print_term_o(t, tout) << std::endl; // TRACE("dioph_eq", tout << "after process_q_with_S\n t:"; print_term_o(t, tout) << std::endl;
tout << "dep:"; print_dep(tout, dep) << std::endl;); // tout << "dep:"; print_dep(tout, dep) << std::endl;);
mpq g = gcd_of_coeffs(t); // mpq g = gcd_of_row(t);
if (g.is_one()) { // if (g.is_one()) {
TRACE("dioph_eq", tout << "g is one" << std::endl;); // TRACE("dioph_eq", tout << "g is one" << std::endl;);
return false; // return false;
} else if (g.is_zero()) { // } else if (g.is_zero()) {
handle_constant_term(t, j, dep); // handle_constant_term(t, j, dep);
if (!m_infeas_explanation.empty()) // if (!m_infeas_explanation.empty())
return true; // return true;
return false; // return false;
} // }
return tighten_bounds_for_term(t, g, j, dep); // return tighten_bounds_for_term(t, g, j, dep);
} }
void handle_constant_term(term_o& t, unsigned j, u_dependency* dep) { void handle_constant_term(term_o& t, unsigned j, u_dependency* dep) {
if (t.c().is_zero()) { if (t.c().is_zero()) {
@ -454,57 +486,57 @@ namespace lp {
// returns true if there is a change // returns true if there is a change
// dep comes from the substitution with S // dep comes from the substitution with S
bool tighten_bounds_for_term(term_o& t, const mpq& g, unsigned j, u_dependency* dep) { bool tighten_bounds_for_term(term_o& t, const mpq& g, unsigned j, u_dependency* dep) {
mpq rs; // mpq rs;
bool is_strict; // bool is_strict;
bool change = false; // bool change = false;
u_dependency *b_dep = nullptr; // u_dependency *b_dep = nullptr;
SASSERT(!g.is_zero()); // SASSERT(!g.is_zero());
if (lra.has_upper_bound(j, b_dep, rs, is_strict)) { // if (lra.has_upper_bound(j, b_dep, rs, is_strict)) {
TRACE("dioph_eq", tout << "tighten upper bound for x" << j << " with rs:" << rs << std::endl;); // TRACE("dioph_eq", tout << "tighten upper bound for x" << j << " with rs:" << rs << std::endl;);
rs = (rs - t.c()) / g; // rs = (rs - t.c()) / g;
TRACE("dioph_eq", tout << "tighten upper bound for x" << j << " with rs:" << rs << std::endl;); // TRACE("dioph_eq", tout << "tighten upper bound for x" << j << " with rs:" << rs << std::endl;);
if (!rs.is_int()) { // if (!rs.is_int()) {
tighten_bound_for_term_for_bound_kind(t, g, j, lra.mk_join(dep, b_dep), rs, true); // tighten_bound_for_term_for_bound_kind(t, g, j, lra.mk_join(dep, b_dep), rs, true);
change = true; // change = true;
} // }
} // }
if (lra.has_lower_bound(j, b_dep, rs, is_strict)) { // if (lra.has_lower_bound(j, b_dep, rs, is_strict)) {
TRACE("dioph_eq", tout << "tighten lower bound for x" << j << " with rs:" << rs << std::endl;); // TRACE("dioph_eq", tout << "tighten lower bound for x" << j << " with rs:" << rs << std::endl;);
rs = (rs - t.c()) / g; // rs = (rs - t.c()) / g;
TRACE("dioph_eq", tout << "tighten lower bound for x" << j << " with rs:" << rs << std::endl;); // TRACE("dioph_eq", tout << "tighten lower bound for x" << j << " with rs:" << rs << std::endl;);
if (!rs.is_int()) { // if (!rs.is_int()) {
tighten_bound_for_term_for_bound_kind(t, g, j, lra.mk_join(b_dep, dep), rs, false); // tighten_bound_for_term_for_bound_kind(t, g, j, lra.mk_join(b_dep, dep), rs, false);
change = true; // change = true;
} // }
} // }
return change; // return change;
} }
void tighten_bound_for_term_for_bound_kind(term_o& t, void tighten_bound_for_term_for_bound_kind(term_o& t,
const mpq& g, unsigned j, u_dependency* dep, const mpq & ub, bool upper) { const mpq& g, unsigned j, u_dependency* dep, const mpq & ub, bool upper) {
// ub = (upper_bound(j) - t.c())/g. // // ub = (upper_bound(j) - t.c())/g.
// we have x[j] = t = g*t_+ t.c() <= upper_bound(j), then // // we have x[j] = t = g*t_+ t.c() <= upper_bound(j), then
// t_ <= floor((upper_bound(j) - t.c())/g) = floor(ub) // // t_ <= floor((upper_bound(j) - t.c())/g) = floor(ub)
// // //
// so xj = g*t_+t.c() <= g*floor(ub) + t.c() is new upper bound // // so xj = g*t_+t.c() <= g*floor(ub) + t.c() is new upper bound
// <= can be replaced with >= for lower bounds, with ceil instead of floor // // <= can be replaced with >= for lower bounds, with ceil instead of floor
mpq bound = g * (upper?floor(ub):ceil(ub))+t.c(); // mpq bound = g * (upper?floor(ub):ceil(ub))+t.c();
TRACE("dioph_eq", tout << "upper:" << upper << std::endl; // TRACE("dioph_eq", tout << "upper:" << upper << std::endl;
tout << "new " << (upper? "upper":"lower" ) << "bound:" << bound << std::endl;); // tout << "new " << (upper? "upper":"lower" ) << "bound:" << bound << std::endl;);
dep = lra.mk_join(dep, upper? lra.get_column_upper_bound_witness(j): lra.get_column_lower_bound_witness(j) ); // dep = lra.mk_join(dep, upper? lra.get_column_upper_bound_witness(j): lra.get_column_lower_bound_witness(j) );
SASSERT(upper && bound <= lra.get_upper_bound(j).x || !upper && bound >= lra.get_lower_bound(j).x); // SASSERT(upper && bound <= lra.get_upper_bound(j).x || !upper && bound >= lra.get_lower_bound(j).x);
lconstraint_kind kind = upper? lconstraint_kind::LE: lconstraint_kind::GE; // lconstraint_kind kind = upper? lconstraint_kind::LE: lconstraint_kind::GE;
lra.update_column_type_and_bound(j, kind, bound, dep); // lra.update_column_type_and_bound(j, kind, bound, dep);
TRACE("dioph_eq", // TRACE("dioph_eq",
tout << "new " << (upper? "upper":"lower" ) << "bound:" << bound << std::endl; // tout << "new " << (upper? "upper":"lower" ) << "bound:" << bound << std::endl;
tout << "bound_dep:\n";print_dep(tout, dep) << std::endl;); // tout << "bound_dep:\n";print_dep(tout, dep) << std::endl;);
} }
public:
lia_move check() { lia_move check() {
init(); init();
while(m_f.size()) { while(m_f.size()) {
@ -519,14 +551,14 @@ namespace lp {
rewrite_eqs(); rewrite_eqs();
} }
TRACE("dioph_eq", print_S(tout);); TRACE("dioph_eq", print_S(tout););
lia_move ret = tighten_with_S(); // lia_move ret = tighten_with_S();
if (ret == lia_move::conflict) { // if (ret == lia_move::conflict) {
lra.settings().stats().m_dio_conflicts++; // lra.settings().stats().m_dio_conflicts++;
return lia_move::conflict; // return lia_move::conflict;
} // }
return lia_move::undef; return lia_move::undef;
} }
private:
std::list<unsigned>::iterator pick_eh() { std::list<unsigned>::iterator pick_eh() {
return m_f.begin(); // TODO: make a smarter joice return m_f.begin(); // TODO: make a smarter joice
} }
@ -561,23 +593,23 @@ namespace lp {
} }
} }
std::tuple<mpq, unsigned, int> find_minimal_abs_coeff(const term_o& eh) const { std::tuple<mpq, unsigned, int> find_minimal_abs_coeff(unsigned row_index) const {
bool first = true, first_fresh = true; bool first = true, first_fresh = true;
mpq ahk, ahk_fresh; mpq ahk, ahk_fresh;
unsigned k, k_fresh; unsigned k, k_fresh;
int k_sign, k_sign_fresh; int k_sign, k_sign_fresh;
mpq t; mpq t;
for (const auto& p : eh) { for (const auto& p : m_e_matrix.m_rows[row_index]) {
t = abs(p.coeff()); t = abs(p.coeff());
if (first_fresh || t < ahk_fresh) { if (first_fresh || t < ahk_fresh) {
ahk_fresh = t; ahk_fresh = t;
k_sign_fresh = p.coeff().is_pos() ? 1 : -1; k_sign_fresh = p.coeff().is_pos() ? 1 : -1;
k_fresh = p.j(); k_fresh = p.var();
first_fresh = false; first_fresh = false;
} else if (first || t < ahk) { } else if (first || t < ahk) {
ahk = t; ahk = t;
k_sign = p.coeff().is_pos() ? 1 : -1; k_sign = p.coeff().is_pos() ? 1 : -1;
k = p.j(); k = p.var();
first = false; first = false;
if (ahk.is_one()) if (ahk.is_one())
break; break;
@ -605,61 +637,105 @@ namespace lp {
TRACE("dioph_eq", tout << "subst_term:"; print_term_o(t, tout) << std::endl;); TRACE("dioph_eq", tout << "subst_term:"; print_term_o(t, tout) << std::endl;);
return t; return t;
} }
void eliminate_var_in_f(eprime_pair& e_pair, unsigned k, int k_sign) {
term_o t = get_term_to_subst(e_pair.m_e, k, k_sign); std::ostream& print_e_row(unsigned i, std::ostream& out) {
substitute_var_on_f(k, k_sign, t, e_pair.m_l, -1) ; // -1 is for the index to ignore return print_term_o(get_term_from_e_matrix(i), out);
}
// j is the variable to eliminate, it appears in row e.m_e_matrix with
// coefficient +-1
void eliminate_var_in_f(eprime_entry& e, unsigned j, int j_sign) {
unsigned piv_row_index = e.m_row_index;
auto &column = m_e_matrix.m_columns[j];
int pivot_col_cell_index = -1;
for (unsigned k = 0; k < column.size(); k++) {
if (column[k].var() == piv_row_index) {
pivot_col_cell_index = k;
break;
}
}
SASSERT(pivot_col_cell_index != -1);
if (pivot_col_cell_index != 0) {
// 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_e_matrix.m_rows[piv_row_index][column[0].offset()].offset() = 0;
m_e_matrix.m_rows[c.var()][c.offset()].offset() = pivot_col_cell_index;
}
unsigned f_rows_in_column =0;
for (const auto& c: column) {
if (m_row_status[c.var()] == row_status::F)
f_rows_in_column ++;
}
TRACE("dioph_eq", tout << "f_rows_in_column:" << f_rows_in_column << std::endl;);
while (column.size() > 1 && f_rows_in_column > 0 ) {
auto & c = column.back();
if (m_row_status[c.var()] != row_status::F)
continue;
f_rows_in_column--;
SASSERT(c.var() != piv_row_index);
mpq coeff = m_e_matrix.get_val(c);
TRACE("dioph_eq", tout << "c_row:" << c.var(); print_e_row(c.var(), tout) << std::endl;);
m_e_matrix.pivot_row_to_row_given_cell_with_sign(piv_row_index, c, j, j_sign);
m_c[c.var()] -= j_sign* coeff*m_c[piv_row_index];
TRACE("dioph_eq", tout << "after pivoting c_row:"; print_e_row(c.var(), tout) << std::endl;);
}
} }
// k is the variable to substitute // k is the variable to substitute
void fresh_var_step(unsigned e_index, unsigned k) { void fresh_var_step(unsigned e_index, unsigned k, const mpq& ahk) {
eprime_pair & e_pair = m_eprime[e_index]; eprime_entry & e = m_eprime[e_index];
// step 7 from the paper unsigned h = e.m_row_index;
auto & eh = e_pair.m_e; // backup the term at h
// xt is the fresh variable m_indexed_work_vector.clear();
unsigned xt = m_last_fresh_x_var++; m_indexed_work_vector.resize(lra.column_count());
TRACE("dioph_eq", tout << "introduce fresh xt:" << "x"<< var_str(xt) << std::endl; auto hrow = m_e_matrix.m_rows[h];
tout << "eh:"; print_term_o(eh,tout) << std::endl;); for (const auto& cell : hrow)
/* Let eh = sum (a_i*x_i) + c m_indexed_work_vector.set_value(cell.coeff(), cell.var());
For each i != k, let a_i = a_qi*ahk + a_ri, and let c = c_q * ahk + c_r while (hrow.size() > 0) {
eh = ahk * (x_k + sum {a_qi*x_i) + c_q | i != k }) + sum {a_ri*x_i | i != k} + c_r = 0 auto & c = hrow.back();
xt = x_k + sum {a_qi*x_i) + c_q | i != k }, it will be xt_term m_e_matrix.remove_element(hrow, c);
Then x_k -> - sum {a_qi*x_i) + c_q | i != k } + xt, it will be subs_k
eh = ahk*xt + ...
*/
term_o xt_term;
term_o k_subs;
// copy it aside
const mpq ahk = eh.get_coeff(k);
mpq r, q = machine_div_rem(eh.c(), ahk, r);
xt_term.add_var(k);
xt_term.c() = q;
xt_term.add_monomial(mpq(-1), xt);
k_subs.add_var(xt);
k_subs.c() = - q;
term_o et; //et is the elimination k from eh, which is an optimization
et.add_monomial(ahk, xt);
et.c() = r;
for (const auto & p: eh) {
if (p.j() == k) continue;
q = machine_div_rem(p.coeff(), ahk, r);
xt_term.add_monomial(q, p.j());
k_subs.add_monomial(-q, p.j());
et.add_monomial(r, p.j());
} }
m_eprime[e_index].m_e = et;
// eprime_pair xt_entry = {xt_term, lar_term()}; // 0 for m_l field // step 7 from the paper
eprime_pair xt_entry = {xt_term, nullptr}; // nullptr for the dependency // xt is the fresh variable
m_eprime.push_back(xt_entry); unsigned xt = m_e_matrix.column_count();
TRACE("dioph_eq", tout << "xt_term:"; print_term_o(xt_term, tout) << std::endl; unsigned fresh_row = m_e_matrix.row_count();
tout << "k_subs:"; print_term_o(k_subs, tout) << std::endl; m_e_matrix.add_row(); // for the fresh variable definition
print_eprime_entry(m_eprime.size()-1, tout);); m_e_matrix.add_column(); // the fresh variable itself
substitute_var_on_f(k, 1, k_subs, xt_entry.m_l, e_index); m_row_status.push_back(row_status::S); // adding a new row to S
// Add a new row for the fresh variable definition
// the term to eliminate the fresh variable /* Let eh = sum(ai*xi) + c. For each i != k, let ai = qi*ahk + ri, and let c = c_q * ahk + c_r.
term_o xt_subs = xt_term.clone(); eh = ahk*(x_k + sum{qi*xi|i != k} + c_q) + sum {ri*xi|i!= k} + c_r.
xt_subs.add_monomial(mpq(1), xt); Then -xt + x_k + sum {qi*x_i)| i != k} + c_q will be the fresh row
TRACE("dioph_eq", tout << "xt_subs: x"<< var_str(xt) << " -> "; print_term_o(xt_subs, tout) << std::endl;); eh = ahk*xt + sum {ri*x_i | i != k} + c_r is the row m_e_matrix[e.m_row_index]
m_sigma.insert(xt, xt_subs); */
mpq q, r;
q = machine_div_rem(m_c[h], ahk, r);
m_c[h] = r;
m_c.push_back(q);
m_e_matrix.add_new_element(h, xt, ahk);
m_e_matrix.add_new_element(fresh_row, xt, -mpq(1));
m_e_matrix.add_new_element(fresh_row, k, mpq(1));
for (unsigned i: m_indexed_work_vector.m_index) {
mpq ai = m_indexed_work_vector[i];
if (i == k) continue;
q = machine_div_rem(ai, ahk, r);
if (!r.is_zero())
m_e_matrix.add_new_element(h, i, r);
if (!q.is_zero())
m_e_matrix.add_new_element(fresh_row, i, q);
}
// add entry to S
m_eprime.push_back({fresh_row, term_o(), nullptr});
this->m_s.push_back(fresh_row);
TRACE("dioph_eq", tout << "changed entry:"; print_eprime_entry(e_index, tout)<< std::endl;
tout << "added to S:\n"; print_eprime_entry(m_eprime.size()-1, tout););
eliminate_var_in_f(m_eprime.back(), k, 1);
} }
std::ostream& print_eprime_entry(unsigned i, std::ostream& out) { std::ostream& print_eprime_entry(unsigned i, std::ostream& out) {
@ -667,16 +743,20 @@ namespace lp {
return print_eprime_entry(m_eprime[i], out); return print_eprime_entry(m_eprime[i], out);
} }
std::ostream& print_eprime_entry(const eprime_pair& e, std::ostream& out) { std::ostream& print_eprime_entry(const eprime_entry& e, std::ostream& out) {
out << "{\n"; out << "{\n";
print_term_o(e.m_e, out << "\tm_e:") << "," << std::endl; out << "\trow:" << e.m_row_index << "," << std::endl;
print_dep(out<< "\tm_l:", e.m_l) << "\n"; print_term_o(get_term_from_e_matrix(e.m_row_index), out << "\trow:");
// print_dep(out<< "\tm_l:", e.m_l) << "\n";
out << "}"<< std::endl; out << "}"<< std::endl;
return out; return out;
} }
// k is the index of the index of the variable with the coefficient +-1 that is being substituted // k is the index of the index of the variable with the coefficient +-1 that is being substituted
void move_entry_from_f_to_s(unsigned k, std::list<unsigned>::iterator it) { void move_entry_from_f_to_s(unsigned k, std::list<unsigned>::iterator it) {
SASSERT(m_row_status[*it] == row_status::F);
m_row_status[*it] = row_status::S;
if (k >= m_k2s.size()) { // k is a fresh variable if (k >= m_k2s.size()) { // k is a fresh variable
m_k2s.resize(k+1, -1 ); m_k2s.resize(k+1, -1 );
} }
@ -691,8 +771,7 @@ namespace lp {
auto eh_it = pick_eh(); auto eh_it = pick_eh();
auto& eprime_entry = m_eprime[*eh_it]; auto& eprime_entry = m_eprime[*eh_it];
TRACE("dioph_eq", print_eprime_entry(*eh_it, tout);); TRACE("dioph_eq", print_eprime_entry(*eh_it, tout););
term_o& eh = eprime_entry.m_e; auto [ahk, k, k_sign] = find_minimal_abs_coeff(eprime_entry.m_row_index);
auto [ahk, k, k_sign] = find_minimal_abs_coeff(eh);
TRACE("dioph_eq", tout << "ahk:" << ahk << ", k:" << k << ", k_sign:" << k_sign << std::endl;); TRACE("dioph_eq", tout << "ahk:" << ahk << ", k:" << k << ", k_sign:" << k_sign << std::endl;);
if (ahk.is_one()) { if (ahk.is_one()) {
eprime_entry.m_e.j() = k; eprime_entry.m_e.j() = k;
@ -700,10 +779,10 @@ namespace lp {
move_entry_from_f_to_s(k, eh_it); move_entry_from_f_to_s(k, eh_it);
eliminate_var_in_f(eprime_entry, k , k_sign); eliminate_var_in_f(eprime_entry, k , k_sign);
} else { } else {
fresh_var_step(*eh_it, k); fresh_var_step(*eh_it, k, ahk);
} }
} }
public:
void explain(explanation& ex) { void explain(explanation& ex) {
if (m_conflict_index == UINT_MAX) { if (m_conflict_index == UINT_MAX) {
SASSERT(!(lra.get_status() == lp_status::FEASIBLE || lra.get_status() == lp_status::OPTIMAL)); SASSERT(!(lra.get_status() == lp_status::FEASIBLE || lra.get_status() == lp_status::OPTIMAL));

1
src/math/lp/static_matrix.cpp
View file

@ -59,6 +59,7 @@ template void static_matrix<mpq, numeric_pair<mpq> >::set(unsigned int, unsigned
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> >::pivot_row_to_row_given_cell_with_sign(unsigned int, column_cell&, unsigned int, 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(vector<lp::row_cell<lp::mpq>, true, unsigned int>&, lp::row_cell<lp::mpq>&);
} }

10
src/math/lp/static_matrix.h
View file

@ -68,7 +68,8 @@ class static_matrix
}; };
std::stack<dim> m_stack; std::stack<dim> m_stack;
public: public:
vector<int> m_vector_of_row_offsets;
vector<int> m_work_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;
@ -110,7 +111,7 @@ public:
static_matrix() = default; static_matrix() = default;
// constructor // constructor
static_matrix(unsigned m, unsigned n): m_vector_of_row_offsets(n, -1) { static_matrix(unsigned m, unsigned n): m_work_vector_of_row_offsets(n, -1) {
init_row_columns(m, n); init_row_columns(m, n);
} }
// constructor that copies columns of the basis from A // constructor that copies columns of the basis from A
@ -133,7 +134,7 @@ public:
void add_row() {m_rows.push_back(row_strip<T>());} void add_row() {m_rows.push_back(row_strip<T>());}
void add_column() { void add_column() {
m_columns.push_back(column_strip()); m_columns.push_back(column_strip());
m_vector_of_row_offsets.push_back(-1); m_work_vector_of_row_offsets.push_back(-1);
} }
void forget_last_columns(unsigned how_many_to_forget); void forget_last_columns(unsigned how_many_to_forget);
@ -289,7 +290,8 @@ public:
} }
// pivot row i to row ii // pivot row i to row ii
bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned); bool pivot_row_to_row_given_cell(unsigned i, column_cell& c, unsigned j);
void pivot_row_to_row_given_cell_with_sign(unsigned piv_row_index, column_cell& c, unsigned j, int j_sign);
void scan_row_ii_to_offset_vector(const row_strip<T> & rvals); void scan_row_ii_to_offset_vector(const row_strip<T> & rvals);
void transpose_rows(unsigned i, unsigned ii) { void transpose_rows(unsigned i, unsigned ii) {

82
src/math/lp/static_matrix_def.h
View file

@ -23,6 +23,7 @@ Revision History:
#include <utility> #include <utility>
#include <set> #include <set>
#include "math/lp/static_matrix.h" #include "math/lp/static_matrix.h"
#include "static_matrix.h"
namespace lp { namespace lp {
// each assignment for this matrix should be issued only once!!! // each assignment for this matrix should be issued only once!!!
@ -31,7 +32,7 @@ inline void addmul(mpq& r, mpq const& a, mpq const& b) { r.addmul(a, b); }
template <typename T, typename X> template <typename T, typename X>
void static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) { void static_matrix<T, X>::init_row_columns(unsigned m, unsigned n) {
lp_assert(m_rows.size() == 0 && m_columns.size() == 0); SASSERT(m_rows.size() == 0 && m_columns.size() == 0);
for (unsigned i = 0; i < m; i++) { for (unsigned i = 0; i < m; i++) {
m_rows.push_back(row_strip<T>()); m_rows.push_back(row_strip<T>());
} }
@ -43,15 +44,16 @@ 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.size(); j++)
m_vector_of_row_offsets[rvals[j].var()] = j; m_work_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.var(); unsigned ii = c.var();
lp_assert(i < row_count() && ii < column_count() && i != ii); SASSERT(i < row_count() && ii < column_count() && i != ii);
T alpha = -get_val(c); T alpha = -get_val(c);
lp_assert(!is_zero(alpha)); SASSERT(!is_zero(alpha));
auto & rowii = m_rows[ii]; auto & rowii = m_rows[ii];
remove_element(rowii, rowii[c.offset()]); remove_element(rowii, rowii[c.offset()]);
scan_row_ii_to_offset_vector(rowii); scan_row_ii_to_offset_vector(rowii);
@ -60,8 +62,8 @@ template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_giv
for (const auto & iv : m_rows[i]) { for (const auto & iv : m_rows[i]) {
unsigned j = iv.var(); unsigned j = iv.var();
if (j == pivot_col) continue; if (j == pivot_col) continue;
lp_assert(!is_zero(iv.coeff())); SASSERT(!is_zero(iv.coeff()));
int j_offs = m_vector_of_row_offsets[j]; int j_offs = m_work_vector_of_row_offsets[j];
if (j_offs == -1) { // it is a new element if (j_offs == -1) { // it is a new element
T alv = alpha * iv.coeff(); T alv = alpha * iv.coeff();
add_new_element(ii, j, alv); add_new_element(ii, j, alv);
@ -72,7 +74,7 @@ template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_giv
} }
// 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].var()] = -1; m_work_vector_of_row_offsets[rowii[k].var()] = -1;
} }
// remove zeroes // remove zeroes
@ -83,11 +85,47 @@ template <typename T, typename X> bool static_matrix<T, X>::pivot_row_to_row_giv
return !rowii.empty(); return !rowii.empty();
} }
template <typename T, typename X>
inline void static_matrix<T, X>::pivot_row_to_row_given_cell_with_sign(unsigned piv_row_index,
column_cell& c, unsigned pivot_col, int pivot_sign) {
unsigned ii = c.var();
SASSERT(ii != piv_row_index);
T alpha = get_val(c) * pivot_sign;
SASSERT(!is_zero(alpha));
auto & rowii = m_rows[ii];
remove_element(rowii, rowii[c.offset()]);
scan_row_ii_to_offset_vector(rowii);
unsigned prev_size_ii = rowii.size();
// run over the pivot row and update row ii
for (const auto & iv : m_rows[piv_row_index]) {
unsigned j = iv.var();
if (j == pivot_col) continue;
SASSERT(!is_zero(iv.coeff()));
int j_offs = m_work_vector_of_row_offsets[j];
if (j_offs == -1) { // it is a new element
T alv = alpha * iv.coeff();
add_new_element(ii, j, alv);
}
else {
addmul(rowii[j_offs].coeff(), iv.coeff(), alpha);
}
}
// clean the work vector
for (unsigned k = 0; k < prev_size_ii; k++) {
m_work_vector_of_row_offsets[rowii[k].var()] = -1;
}
// remove zeroes
for (unsigned k = rowii.size(); k-- > 0; ) {
if (is_zero(rowii[k].coeff()))
remove_element(rowii, rowii[k]);
}
}
// constructor that copies columns of the basis from A // constructor that copies columns of the basis from A
template <typename T, typename X> template <typename T, typename X>
static_matrix<T, X>::static_matrix(static_matrix const &A, unsigned * /* basis */) : static_matrix<T, X>::static_matrix(static_matrix const &A, unsigned * /* basis */) :
m_vector_of_row_offsets(A.column_count(), numeric_traits<T>::zero()) { m_work_vector_of_row_offsets(A.column_count(), numeric_traits<T>::zero()) {
unsigned m = A.row_count(); unsigned m = A.row_count();
init_row_columns(m, m); init_row_columns(m, m);
for (; m-- > 0; ) for (; m-- > 0; )
@ -96,14 +134,14 @@ static_matrix<T, X>::static_matrix(static_matrix const &A, unsigned * /* basis *
} }
template <typename T, typename X> void static_matrix<T, X>::clear() { template <typename T, typename X> void static_matrix<T, X>::clear() {
m_vector_of_row_offsets.clear(); m_work_vector_of_row_offsets.clear();
m_rows.clear(); m_rows.clear();
m_columns.clear(); m_columns.clear();
} }
template <typename T, typename X> void static_matrix<T, X>::init_vector_of_row_offsets() { template <typename T, typename X> void static_matrix<T, X>::init_vector_of_row_offsets() {
m_vector_of_row_offsets.clear(); m_work_vector_of_row_offsets.clear();
m_vector_of_row_offsets.resize(column_count(), -1); m_work_vector_of_row_offsets.resize(column_count(), -1);
} }
template <typename T, typename X> void static_matrix<T, X>::init_empty_matrix(unsigned m, unsigned n) { template <typename T, typename X> void static_matrix<T, X>::init_empty_matrix(unsigned m, unsigned n) {
@ -112,9 +150,9 @@ template <typename T, typename X> void static_matrix<T, X>::init_empty_matrix(un
} }
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()); SASSERT(col < column_count());
column_strip & colstrip = m_columns[col]; column_strip & colstrip = m_columns[col];
lp_assert(colstrip.size() > 0); SASSERT(colstrip.size() > 0);
unsigned ret = 0; unsigned ret = 0;
for (auto & t : colstrip) { for (auto & t : colstrip) {
if (t.var() > ret) { if (t.var() > ret) {
@ -125,7 +163,7 @@ template <typename T, typename X> unsigned static_matrix<T, X>::lowest_row_in_co
} }
template <typename T, typename X> void static_matrix<T, X>::forget_last_columns(unsigned how_many_to_forget) { template <typename T, typename X> void static_matrix<T, X>::forget_last_columns(unsigned how_many_to_forget) {
lp_assert(m_columns.size() >= how_many_to_forget); SASSERT(m_columns.size() >= how_many_to_forget);
unsigned j = column_count() - 1; unsigned j = column_count() - 1;
for (; how_many_to_forget-- > 0; ) { for (; how_many_to_forget-- > 0; ) {
remove_last_column(j --); remove_last_column(j --);
@ -146,12 +184,12 @@ template <typename T, typename X> void static_matrix<T, X>::remove_last_column(u
} }
} }
m_columns.pop_back(); m_columns.pop_back();
m_vector_of_row_offsets.pop_back(); m_work_vector_of_row_offsets.pop_back();
} }
template <typename T, typename X> void static_matrix<T, X>::set(unsigned row, unsigned col, T const & val) { 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; if (numeric_traits<T>::is_zero(val)) return;
lp_assert(row < row_count() && col < column_count()); SASSERT(row < row_count() && col < column_count());
auto & r = m_rows[row]; auto & r = m_rows[row];
unsigned offs_in_cols = m_columns[col].size(); unsigned offs_in_cols = m_columns[col].size();
m_columns[col].push_back(make_column_cell(row, r.size())); m_columns[col].push_back(make_column_cell(row, r.size()));
@ -229,7 +267,7 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency()
for (unsigned i = 0; i < m_rows.size(); i++) { for (unsigned 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.var()); std::pair<unsigned, unsigned> p(i, t.var());
lp_assert(by_rows.find(p) == by_rows.end()); SASSERT(by_rows.find(p) == by_rows.end());
by_rows[p] = t.coeff(); by_rows[p] = t.coeff();
} }
} }
@ -237,16 +275,16 @@ template <typename T, typename X> void static_matrix<T, X>::check_consistency()
for (unsigned i = 0; i < m_columns.size(); i++) { for (unsigned 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.var(), i); std::pair<unsigned, unsigned> p(t.var(), i);
lp_assert(by_cols.find(p) == by_cols.end()); SASSERT(by_cols.find(p) == by_cols.end());
by_cols[p] = get_val(t); by_cols[p] = get_val(t);
} }
} }
lp_assert(by_rows.size() == by_cols.size()); SASSERT(by_rows.size() == by_cols.size());
for (auto & t : by_rows) { for (auto & t : by_rows) {
auto ic = by_cols.find(t.first); auto ic = by_cols.find(t.first);
lp_assert(ic != by_cols.end()); SASSERT(ic != by_cols.end());
lp_assert(t.second == ic->second); SASSERT(t.second == ic->second);
} }
} }
#endif #endif