From ad2445e42337c7ab347a3e3b207603cd12055473 Mon Sep 17 00:00:00 2001 From: Nikolaj Bjorner Date: Mon, 9 May 2022 16:33:05 -0700 Subject: [PATCH] gauss jordan Signed-off-by: Nikolaj Bjorner --- src/math/simplex/simplex.cpp | 4 ++ src/math/simplex/simplex.h | 2 + src/math/simplex/sparse_matrix.h | 69 ++++++++++++++++++++++++++-- src/math/simplex/sparse_matrix_ops.h | 68 +++++++++++++++++++++++++++ src/sat/smt/arith_diagnostics.cpp | 13 ++++++ src/sat/smt/arith_solver.cpp | 2 +- src/sat/smt/arith_solver.h | 3 ++ src/test/simplex.cpp | 31 +++++++++++-- 8 files changed, 183 insertions(+), 9 deletions(-) create mode 100644 src/math/simplex/sparse_matrix_ops.h diff --git a/src/math/simplex/simplex.cpp b/src/math/simplex/simplex.cpp index 6d1c4d348..7649586df 100644 --- a/src/math/simplex/simplex.cpp +++ b/src/math/simplex/simplex.cpp @@ -18,6 +18,7 @@ Notes: --*/ #include "math/simplex/simplex.h" +#include "math/simplex/sparse_matrix_ops.h" #include "math/simplex/sparse_matrix_def.h" #include "math/simplex/simplex_def.h" #include "util/rational.h" @@ -36,6 +37,9 @@ namespace simplex { } } + void gauss_jordan(sparse_matrix& M) { + sparse_matrix_ops::gauss_jordan(M); + } void ensure_rational_solution(simplex& S) { rational delta(1); diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h index b8781d4ab..44ff3191e 100644 --- a/src/math/simplex/simplex.h +++ b/src/math/simplex/simplex.h @@ -200,5 +200,7 @@ namespace simplex { }; void ensure_rational_solution(simplex& s); + + void gauss_jordan(sparse_matrix& s); }; diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h index f18bb6f87..ffbe009e0 100644 --- a/src/math/simplex/sparse_matrix.h +++ b/src/math/simplex/sparse_matrix.h @@ -201,8 +201,22 @@ namespace simplex { row_iterator row_begin(row const& r) { return row_iterator(m_rows[r.id()], true); } row_iterator row_end(row const& r) { return row_iterator(m_rows[r.id()], false); } + class row_vars { + friend class sparse_matrix; + sparse_matrix& s; + row r; + row_vars(sparse_matrix& s, row r): s(s), r(r) {} + public: + row_iterator begin() { return s.row_begin(r); } + row_iterator end() { return s.row_end(r); } + }; + + row_vars get_row(row r) { return row_vars(*this, r); } + unsigned column_size(var_t v) const { return m_columns[v].size(); } + unsigned num_vars() const { return m_columns.size(); } + class col_iterator { friend class sparse_matrix; unsigned m_curr; @@ -228,15 +242,16 @@ namespace simplex { --m_col.m_refs; } - row get_row() { + row get_row() const { return row(m_col.m_entries[m_curr].m_row_id); } - row_entry const& get_row_entry() { + row_entry const& get_row_entry() const { col_entry const& c = m_col.m_entries[m_curr]; int row_id = c.m_row_id; return m_rows[row_id].m_entries[c.m_row_idx]; } - + + std::pair operator*() { return std::make_pair(get_row(), &get_row_entry()); } col_iterator & operator++() { ++m_curr; move_to_used(); return *this; } col_iterator operator++(int) { col_iterator tmp = *this; ++*this; return tmp; } bool operator==(col_iterator const & it) const { return m_curr == it.m_curr; } @@ -246,10 +261,58 @@ namespace simplex { col_iterator col_begin(int v) const { return col_iterator(m_columns[v], m_rows, true); } col_iterator col_end(int v) const { return col_iterator(m_columns[v], m_rows, false); } + class var_rows { + friend class sparse_matrix; + sparse_matrix& s; + int v; + var_rows(sparse_matrix& s, int v):s(s), v(v) {} + public: + col_iterator begin() { return s.col_begin(v); } + col_iterator end() { return s.col_end(v); } + }; + + var_rows get_rows(int v) { return var_rows(*this, v); } + + class all_row_iterator { + friend class sparse_matrix; + unsigned m_curr; + vector<_row> const& m_rows; + void move_to_next() { + while (m_curr < m_rows.size() && m_rows[m_curr].size() == 0) { + std::cout << "size is 0 for " << m_curr << "\n"; + ++m_curr; + } + } + public: + all_row_iterator(unsigned curr, vector<_row> const& rows): m_curr(curr), m_rows(rows) { + move_to_next(); + } + row operator*() { return row(m_curr); } + all_row_iterator & operator++() { m_curr++; move_to_next(); return *this; } + all_row_iterator operator++(int) { all_row_iterator tmp = *this; ++*this; return tmp; } + bool operator==(all_row_iterator const& it) const { return m_curr == it.m_curr; } + bool operator!=(all_row_iterator const& it) const { return m_curr != it.m_curr; } + }; + + class all_rows { + friend class sparse_matrix; + sparse_matrix& s; + all_rows(sparse_matrix& s): s(s) {} + public: + all_row_iterator begin() { return all_row_iterator(0, s.m_rows); } + all_row_iterator end() { return all_row_iterator(s.m_rows.size(), s.m_rows); } + }; + + + all_rows get_rows() { return all_rows(*this); } + + void display(std::ostream& out); void display_row(std::ostream& out, row const& r); bool well_formed() const; + manager& get_manager() { return m; } + void collect_statistics(::statistics & st) const; }; diff --git a/src/math/simplex/sparse_matrix_ops.h b/src/math/simplex/sparse_matrix_ops.h new file mode 100644 index 000000000..249c9bc1e --- /dev/null +++ b/src/math/simplex/sparse_matrix_ops.h @@ -0,0 +1,68 @@ +/*++ +Copyright (c) 2014 Microsoft Corporation + +Module Name: + + sparse_matrix_ops.h + +Abstract: + + +Author: + + Nikolaj Bjorner (nbjorner) 2014-01-15 + +Notes: + +--*/ + +#pragma once + +#include "math/simplex/sparse_matrix.h" +#include "util/rational.h" + +namespace simplex { + + class sparse_matrix_ops { + public: + static void gauss_jordan(sparse_matrix& M) { + mpq_ext::numeral coeff; + vector c, d; + unsigned m = M.num_vars(); + for (unsigned v = 0; v < m; ++v) + c.push_back(0); + + for (auto const& row : M.get_rows()) { + // scan for non-zero variable in row + bool found = false; + for (auto const& [coeff1, v] : M.get_row(row)) { + if (mpq_manager::is_zero(coeff1)) + continue; + found = true; + d.push_back(v); + c[v] = row.id() + 1; + // eliminate v from other rows. + for (auto const& [row2, row_entry2] : M.get_rows(v)) { + if (row.id() == row2.id()) + continue; + if (row_entry2->m_coeff == 0) + continue; + M.get_manager().set(coeff, (- row_entry2->m_coeff / coeff1).to_mpq()); + M.add(row2, coeff, row); + } + break; + } + if (!found) + d.push_back(0); + + } + + M.get_manager().del(coeff); + + // TODO: do something with c and d + } + }; + + +} + diff --git a/src/sat/smt/arith_diagnostics.cpp b/src/sat/smt/arith_diagnostics.cpp index 2dfd508d1..28105472d 100644 --- a/src/sat/smt/arith_diagnostics.cpp +++ b/src/sat/smt/arith_diagnostics.cpp @@ -79,4 +79,17 @@ namespace arith { lp().settings().stats().collect_statistics(st); if (m_nla) m_nla->collect_statistics(st); } + + char const* solver::bounds_pragma() { + if (!ctx.use_drat()) + return nullptr; + m_bounds_pragma.clear(); + m_bounds_pragma += "bounds "; + for (sat::literal c : m_core) { + if (c.sign()) m_bounds_pragma += "-"; + m_bounds_pragma += std::to_string(c.var()); + m_bounds_pragma += " "; + } + return m_bounds_pragma.c_str(); + } } diff --git a/src/sat/smt/arith_solver.cpp b/src/sat/smt/arith_solver.cpp index 151bc81d1..14ed76738 100644 --- a/src/sat/smt/arith_solver.cpp +++ b/src/sat/smt/arith_solver.cpp @@ -255,7 +255,7 @@ namespace arith { TRACE("arith", for (auto lit : m_core) tout << lit << ": " << s().value(lit) << "\n";); DEBUG_CODE(for (auto lit : m_core) { VERIFY(s().value(lit) == l_true); }); ++m_stats.m_bound_propagations1; - assign(lit, m_core, m_eqs, "bounds"); + assign(lit, m_core, m_eqs, bounds_pragma()); } if (should_refine_bounds() && first) diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h index 3496e42d3..5906ed736 100644 --- a/src/sat/smt/arith_solver.h +++ b/src/sat/smt/arith_solver.h @@ -419,6 +419,9 @@ namespace arith { void false_case_of_check_nla(const nla::lemma& l); void dbg_finalize_model(model& mdl); + std::string m_bounds_pragma; + char const* bounds_pragma(); + public: solver(euf::solver& ctx, theory_id id); diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp index 713cdd31c..e66a880d6 100644 --- a/src/test/simplex.cpp +++ b/src/test/simplex.cpp @@ -15,6 +15,7 @@ Copyright (c) 2015 Microsoft Corporation #define R rational typedef simplex::simplex Simplex; typedef simplex::sparse_matrix sparse_matrix; +typedef simplex::sparse_matrix qmatrix; static vector vec(int i, int j) { vector nv; @@ -24,11 +25,11 @@ static vector vec(int i, int j) { return nv; } -// static vector vec(int i, int j, int k) { -// vector nv = vec(i, j); -// nv.push_back(R(k)); -// return nv; -// } +static vector vec(int i, int j, int k) { + vector nv = vec(i, j); + nv.push_back(R(k)); + return nv; +} // static vector vec(int i, int j, int k, int l) { // vector nv = vec(i, j, k); @@ -131,6 +132,25 @@ static void test4() { feas(S); } +static void add(qmatrix& m, vector const& v) { + m.ensure_var(v.size()); + auto r = m.mk_row(); + for (unsigned u = 0; u < v.size(); ++u) + m.add_var(r, v[u].to_mpq(), u); +} + +static void test5() { + unsynch_mpq_manager m; + qmatrix M(m); + add(M, vec(1, 2, 3)); + add(M, vec(2, 2, 4)); + M.display(std::cout); + gauss_jordan(M); + std::cout << "after\n"; + M.display(std::cout); + +} + void tst_simplex() { reslimit rl; Simplex S(rl); @@ -166,4 +186,5 @@ void tst_simplex() { test2(); test3(); test4(); + test5(); }