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 void sparse_matrix::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 void sparse_matrix::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& M, vector>& K) { - rational D; - vector 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 + static void kernel(sparse_matrix &M, vector> &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::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 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::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()); - for (unsigned i = 0; i < n; ++i) { - if (d[i] > 0) { - auto r = sparse_matrix::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()); + for (unsigned i = 0; i < n_vars; ++i) { + if (d[i] > 0) { + auto r = sparse_matrix::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 &M, vector> &K) { + kernel(M, K); + } +}; +} // namespace simplex