gauss jordan

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

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<mpq_ext>& M) {
+        sparse_matrix_ops::gauss_jordan(M);
+    }
 
     void ensure_rational_solution(simplex<mpq_ext>& S) {
         rational delta(1);

src/math/simplex/simplex.h

@@ -200,5 +200,7 @@ namespace simplex {
     };
 
     void ensure_rational_solution(simplex<mpq_ext>& s);
+
+    void gauss_jordan(sparse_matrix<mpq_ext>& s);
 };
 

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<row, row_entry const*> 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;
 
     };

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<mpq_ext>& M) {
+            mpq_ext::numeral coeff;
+            vector<unsigned> 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<false>::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
+        }
+    };
+
+
+}
+

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();
+    }
 }

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)

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);

src/test/simplex.cpp

@@ -15,6 +15,7 @@ Copyright (c) 2015 Microsoft Corporation
 #define R rational
 typedef simplex::simplex<simplex::mpz_ext> Simplex;
 typedef simplex::sparse_matrix<simplex::mpz_ext> sparse_matrix;
+typedef simplex::sparse_matrix<simplex::mpq_ext> qmatrix;
 
 static vector<R> vec(int i, int j) {
     vector<R> nv;
@@ -24,11 +25,11 @@ static vector<R> vec(int i, int j) {
     return nv;
 }
 
-// static vector<R> vec(int i, int j, int k) {
-//     vector<R> nv = vec(i, j);
-//     nv.push_back(R(k));
-//     return nv;
-// }
+static vector<R> vec(int i, int j, int k) {
+     vector<R> nv = vec(i, j);
+     nv.push_back(R(k));
+     return nv;
+}
 
 // static vector<R> vec(int i, int j, int k, int l) {
 //     vector<R> nv = vec(i, j, k);
@@ -131,6 +132,25 @@ static void test4() {
     feas(S);
 }
 
+static void add(qmatrix& m, vector<R> 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();
 }