adding K

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

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

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

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

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