From 7f62fa2b669864e4107972ab2b619da8b87e175b Mon Sep 17 00:00:00 2001
From: Arie Gurfinkel <arie.gurfinkel@gmail.com>
Date: Fri, 13 May 2022 17:30:35 -0700
Subject: [PATCH] Sparse matrix kernel (#6035)

* Subtle bug in kernel computation

Coefficient was being passed by reference and, therefore, was
being changed indirectly.

In the process, updated the code to be more generic to avoid rational
computation in the middle of matrix manipulation.

* sparse_matrix: fixed handling of 0 in add_var() and add()

particularly in add_var(), without the fix the user is responsible for checking
coefficients for 0.
---
 src/math/simplex/sparse_matrix_def.h |   2 +
 src/math/simplex/sparse_matrix_ops.h | 106 ++++++++++++++-------------
 2 files changed, 57 insertions(+), 51 deletions(-)

diff --git a/src/math/simplex/sparse_matrix_def.h b/src/math/simplex/sparse_matrix_def.h
index a3170461b..b20aa76ba 100644
--- a/src/math/simplex/sparse_matrix_def.h
+++ b/src/math/simplex/sparse_matrix_def.h
@@ -299,6 +299,7 @@ namespace simplex {
 
     template<typename Ext>
     void sparse_matrix<Ext>::add_var(row dst, numeral const& n, var_t v) {
+        if (m.is_zero(n))  return;
         _row& r   = m_rows[dst.id()];
         column& c = m_columns[v];
         unsigned r_idx;
@@ -317,6 +318,7 @@ namespace simplex {
     */
     template<typename Ext>
     void sparse_matrix<Ext>::add(row row1, numeral const& n, row row2) {
+        if (m.is_zero(n)) return;
         m_stats.m_add_rows++;
         _row & r1 = m_rows[row1.id()];
         
diff --git a/src/math/simplex/sparse_matrix_ops.h b/src/math/simplex/sparse_matrix_ops.h
index 8ed36f8bb..039d2fede 100644
--- a/src/math/simplex/sparse_matrix_ops.h
+++ b/src/math/simplex/sparse_matrix_ops.h
@@ -7,7 +7,7 @@ Module Name:
 
 Abstract:
 
-    
+
 Author:
 
     Nikolaj Bjorner (nbjorner) 2014-01-15
@@ -23,59 +23,63 @@ Notes:
 
 namespace simplex {
 
-    class sparse_matrix_ops {
-    public:
-        static void kernel(sparse_matrix<mpq_ext>& M, vector<vector<rational>>& K) {
-            rational D;
-            vector<unsigned> d, c;
-            unsigned n = M.num_vars(), m = M.num_rows();
-            auto& mgr = M.get_manager();
-            c.resize(m, 0u);
-            d.resize(n, 0u);
+class sparse_matrix_ops {
+  public:
+    template <typename Ext>
+    static void kernel(sparse_matrix<Ext> &M, vector<vector<rational>> &K) {
+        using scoped_numeral = typename Ext::scoped_numeral;
 
-            for (unsigned k = 0; k < n; ++k) {
-                d[k] = 0;
-                for (auto [row, row_entry] : M.get_rows(k)) {
-                    if (c[row.id()] != 0)
-                        continue;                    
-                    auto& m_jk = row_entry->m_coeff;
-                    if (mpq_manager<false>::is_zero(m_jk))
-                        continue;
-                    D = rational(-1) / m_jk;
-                    M.mul(row, D.to_mpq());
-                    for (auto [row_i, row_i_entry] : M.get_rows(k)) {
-                        if (row_i.id() == row.id())
-                            continue;                        
-                        auto& m_ik = row_i_entry->m_coeff;
-                        // row_i += m_ik * row
-                        M.add(row_i, m_ik, row);                        
-                    }                    
-                    c[row.id()] = k + 1;
-                    d[k] = row.id() + 1;
-                    break;
+        vector<unsigned> d, c;
+        unsigned n_vars = M.num_vars(), n_rows = M.num_rows();
+        c.resize(n_rows, 0u);
+        d.resize(n_vars, 0u);
+
+        auto &m = M.get_manager();
+        scoped_numeral m_ik(m);
+        scoped_numeral D(m);
+
+        for (unsigned k = 0; k < n_vars; ++k) {
+            d[k] = 0;
+            for (auto [row, row_entry] : M.get_rows(k)) {
+                if (c[row.id()] != 0) continue;
+                auto &m_jk = row_entry->m_coeff;
+                if (mpq_manager<false>::is_zero(m_jk)) continue;
+
+                // D = rational(-1) / m_jk;
+                m.set(D, m_jk);
+                m.inv(D);
+                m.neg(D);
+
+                M.mul(row, D);
+                for (auto [row_i, row_i_entry] : M.get_rows(k)) {
+                    if (row_i.id() == row.id()) continue;
+                    m.set(m_ik, row_i_entry->m_coeff);
+                    // row_i += m_ik * row
+                    M.add(row_i, m_ik, row);
                 }
+                c[row.id()] = k + 1;
+                d[k] = row.id() + 1;
+                break;
             }
-
-            for (unsigned k = 0; k < n; ++k) {
-                if (d[k] != 0)
-                    continue;
-                K.push_back(vector<rational>());
-                for (unsigned i = 0; i < n; ++i) {
-                    if (d[i] > 0) {
-                        auto r = sparse_matrix<mpq_ext>::row(d[i]-1);
-                        K.back().push_back(rational(M.get_coeff(r, k)));
-                    }
-                    else if (i == k)
-                        K.back().push_back(rational(1));
-                    else
-                        K.back().push_back(rational(0));
-                }
-            }
-
         }
 
-    };
-
-
-}
+        for (unsigned k = 0; k < n_vars; ++k) {
+            if (d[k] != 0) continue;
+            K.push_back(vector<rational>());
+            for (unsigned i = 0; i < n_vars; ++i) {
+                if (d[i] > 0) {
+                    auto r = sparse_matrix<mpq_ext>::row(d[i] - 1);
+                    K.back().push_back(rational(M.get_coeff(r, k)));
+                } else if (i == k)
+                    K.back().push_back(rational(1));
+                else
+                    K.back().push_back(rational(0));
+            }
+        }
+    }
 
+    static void kernel(sparse_matrix<mpq_ext> &M, vector<vector<rational>> &K) {
+        kernel<mpq_ext>(M, K);
+    }
+};
+} // namespace simplex