From 6a8ac5f9b195cd5d0ad92608c598848d2c6e7ff0 Mon Sep 17 00:00:00 2001
From: Nikolaj Bjorner <nbjorner@microsoft.com>
Date: Mon, 9 May 2022 16:47:26 -0700
Subject: [PATCH] adding K

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

diff --git a/src/math/simplex/simplex.cpp b/src/math/simplex/simplex.cpp
index 7649586df..372bbc221 100644
--- a/src/math/simplex/simplex.cpp
+++ b/src/math/simplex/simplex.cpp
@@ -37,8 +37,8 @@ namespace simplex {
         }
     }
 
-    void gauss_jordan(sparse_matrix<mpq_ext>& M) {
-        sparse_matrix_ops::gauss_jordan(M);
+    void kernel(sparse_matrix<mpq_ext>& M, vector<vector<mpq>>& K) {
+        sparse_matrix_ops::kernel(M, K);
     }
 
     void ensure_rational_solution(simplex<mpq_ext>& S) {
diff --git a/src/math/simplex/simplex.h b/src/math/simplex/simplex.h
index 44ff3191e..58efef553 100644
--- a/src/math/simplex/simplex.h
+++ b/src/math/simplex/simplex.h
@@ -201,6 +201,6 @@ namespace simplex {
 
     void ensure_rational_solution(simplex<mpq_ext>& s);
 
-    void gauss_jordan(sparse_matrix<mpq_ext>& s);
+    void kernel(sparse_matrix<mpq_ext>& s, vector<vector<mpq>>& K);
 };
 
diff --git a/src/math/simplex/sparse_matrix_ops.h b/src/math/simplex/sparse_matrix_ops.h
index 249c9bc1e..5545f5774 100644
--- a/src/math/simplex/sparse_matrix_ops.h
+++ b/src/math/simplex/sparse_matrix_ops.h
@@ -25,22 +25,25 @@ namespace simplex {
 
     class sparse_matrix_ops {
     public:
-        static void gauss_jordan(sparse_matrix<mpq_ext>& M) {
+        static void kernel(sparse_matrix<mpq_ext>& M, vector<vector<mpq>>& K) {
             mpq_ext::numeral coeff;
-            vector<unsigned> c, d;
+            vector<unsigned> d;
+            bool_vector c;
             unsigned m = M.num_vars();
             for (unsigned v = 0; v < m; ++v)
-                c.push_back(0);
+                c.push_back(false);
             
             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)) {
                     if (mpq_manager<false>::is_zero(coeff1))
                         continue;
-                    found = true;
-                    d.push_back(v);
-                    c[v] = row.id() + 1;
+                    if (c[v])
+                        continue;
+                    d.back() = v + 1;
+                    c[v] = true;
                     // eliminate v from other rows.
                     for (auto const& [row2, row_entry2] : M.get_rows(v)) {
                         if (row.id() == row2.id())
@@ -51,15 +54,21 @@ namespace simplex {
                         M.add(row2, coeff, row);
                     }
                     break;
-                }
-                if (!found)
-                    d.push_back(0);
-                
+                }                
             }
 
             M.get_manager().del(coeff);
 
-            // TODO: do something with c and d
+            // TODO: extract kernel using d
+            for (unsigned k = 0; k < d.size(); ++k) {
+                if (d[k] != 0)
+                    continue;
+                K.push_back(vector<mpq>());
+                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 e66a880d6..f351bc15e 100644
--- a/src/test/simplex.cpp
+++ b/src/test/simplex.cpp
@@ -145,7 +145,8 @@ static void test5() {
     add(M, vec(1, 2, 3));
     add(M, vec(2, 2, 4));    
     M.display(std::cout);
-    gauss_jordan(M);
+    vector<vector<mpq>> K;
+    kernel(M, K);
     std::cout << "after\n";
     M.display(std::cout);