From 0557d72d1c5e7d3ea42eae41dddcdce58cca3627 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Tue, 10 May 2022 07:42:32 -0700
Subject: [PATCH] na

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>
---
 src/math/simplex/simplex.cpp         |  2 +-
 src/math/simplex/simplex.h           |  3 +-
 src/math/simplex/sparse_matrix_ops.h | 41 +++++++++++++++++++---------
 src/test/simplex.cpp                 |  2 +-
 4 files changed, 32 insertions(+), 16 deletions(-)

diff --git a/src/math/simplex/simplex.cpp b/src/math/simplex/simplex.cpp
index 372bbc221..643db9c5b 100644
--- a/src/math/simplex/simplex.cpp
+++ b/src/math/simplex/simplex.cpp
@@ -37,7 +37,7 @@ namespace simplex {
         }
     }
 
-    void kernel(sparse_matrix<mpq_ext>& M, vector<vector<mpq>>& K) {
+    void kernel(sparse_matrix<mpq_ext>& M, vector<vector<rational>>& K) {
         sparse_matrix_ops::kernel(M, K);
     }
 
diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h
index 58efef553..8c470a242 100644
--- a/src/math/simplex/simplex.h
+++ b/src/math/simplex/simplex.h
@@ -33,6 +33,7 @@ Notes:
 
 #include "math/simplex/sparse_matrix.h"
 #include "util/mpq_inf.h"
+#include "util/rational.h"
 #include "util/heap.h"
 #include "util/lbool.h"
 #include "util/uint_set.h"
@@ -201,6 +202,6 @@ namespace simplex {
 
     void ensure_rational_solution(simplex<mpq_ext>& s);
 
-    void kernel(sparse_matrix<mpq_ext>& s, vector<vector<mpq>>& K);
+    void kernel(sparse_matrix<mpq_ext>& s, vector<vector<rational>>& K);
 };
 
diff --git a/src/math/simplex/sparse_matrix_ops.h b/src/math/simplex/sparse_matrix_ops.h
index 5545f5774..e12f421a1 100644
--- a/src/math/simplex/sparse_matrix_ops.h
+++ b/src/math/simplex/sparse_matrix_ops.h
@@ -25,45 +25,60 @@ namespace simplex {
 
     class sparse_matrix_ops {
     public:
-        static void kernel(sparse_matrix<mpq_ext>& M, vector<vector<mpq>>& K) {
+        static void kernel(sparse_matrix<mpq_ext>& M, vector<vector<rational>>& K) {
             mpq_ext::numeral coeff;
-            vector<unsigned> d;
-            bool_vector c;
+            rational D1, D2;
+            vector<unsigned> d, c;
             unsigned m = M.num_vars();
+            auto& mgr = M.get_manager();
             for (unsigned v = 0; v < m; ++v)
-                c.push_back(false);
+                c.push_back(0);
             
             for (auto const& row : M.get_rows()) {
                 // scan for non-zero variable in row
                 bool found = false;
                 d.push_back(0);
-                for (auto const& [coeff1, v] : M.get_row(row)) {
+                for (auto& [coeff1, v] : M.get_row(row)) {
                     if (mpq_manager<false>::is_zero(coeff1))
                         continue;
-                    if (c[v])
+                    if (c[v] != 0)
                         continue;
                     d.back() = v + 1;
-                    c[v] = true;
+                    c[v] = row.id() + 1;
+                    D1 = rational(-1) / coeff1;
+                    mgr.set(coeff1, mpq(-1));
                     // eliminate v from other rows.
                     for (auto const& [row2, row_entry2] : M.get_rows(v)) {
-                        if (row.id() == row2.id())
+                        if (row.id() >= row2.id() || row_entry2->m_coeff == 0)
                             continue;
-                        if (row_entry2->m_coeff == 0)
+                        for (auto& [coeff2, w] : M.get_row(row2)) {
+                            if (v == w)
+                                mgr.set(coeff2, (D1*coeff2).to_mpq());
+                        }                        
+                    }
+                    
+                    for (auto& [coeff2, w] : M.get_row(row)) {
+                        if (v == w)
                             continue;
-                        M.get_manager().set(coeff, (- row_entry2->m_coeff / coeff1).to_mpq());
-                        M.add(row2, coeff, row);
+                        D2 = coeff2;
+                        mgr.set(coeff2, mpq(0));
+                        for (auto const& [row2, row_entry2] : M.get_rows(w)) {
+                            if (row.id() >= row2.id() || row_entry2->m_coeff == 0 || row_entry2->m_var == v)
+                                continue;
+                            // mgr.set(row_entry2->m_coeff, row_entry2->m_coeff + D2*row2[v]);
+                        }
                     }
                     break;
                 }                
             }
 
-            M.get_manager().del(coeff);
+            mgr.del(coeff);
 
             // TODO: extract kernel using d
             for (unsigned k = 0; k < d.size(); ++k) {
                 if (d[k] != 0)
                     continue;
-                K.push_back(vector<mpq>());
+                K.push_back(vector<rational>());
                 for (unsigned i = 0; i < d.size(); ++i) {                    
                     // K.back().push_back(d[i] > 0 ? M[d[i]-1][k] : (i == k) ? 1 : 0);
                 }
diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp
index f351bc15e..a3be5d3bb 100644
--- a/src/test/simplex.cpp
+++ b/src/test/simplex.cpp
@@ -145,7 +145,7 @@ static void test5() {
     add(M, vec(1, 2, 3));
     add(M, vec(2, 2, 4));    
     M.display(std::cout);
-    vector<vector<mpq>> K;
+    vector<vector<rational>> K;
     kernel(M, K);
     std::cout << "after\n";
     M.display(std::cout);