diff --git a/src/math/simplex/sparse_matrix.h b/src/math/simplex/sparse_matrix.h index 669b8c2b9..b51b57181 100644 --- a/src/math/simplex/sparse_matrix.h +++ b/src/math/simplex/sparse_matrix.h @@ -136,6 +136,7 @@ namespace simplex { svector m_var_pos; // temporary map from variables to positions in row unsigned_vector m_var_pos_idx; // indices in m_var_pos stats m_stats; + scoped_numeral m_zero; bool well_formed_row(unsigned row_id) const; bool well_formed_column(unsigned column_id) const; @@ -144,7 +145,7 @@ namespace simplex { public: - sparse_matrix(manager& _m): m(_m) {} + sparse_matrix(manager& _m): m(_m), m_zero(m) {} ~sparse_matrix(); void reset(); @@ -216,6 +217,7 @@ namespace simplex { unsigned column_size(var_t v) const { return m_columns[v].size(); } unsigned num_vars() const { return m_columns.size(); } + unsigned num_rows() const { return m_rows.size(); } class col_iterator { friend class sparse_matrix; @@ -306,11 +308,11 @@ namespace simplex { all_rows get_rows() { return all_rows(*this); } - numeral& get_coeff(row r, unsigned v) { + numeral const& get_coeff(row r, unsigned v) { for (auto & [coeff, u] : get_row(r)) if (u == v) return coeff; - throw default_exception("variable not in row"); + return m_zero; } diff --git a/src/math/simplex/sparse_matrix_ops.h b/src/math/simplex/sparse_matrix_ops.h index 5b508afcf..8ed36f8bb 100644 --- a/src/math/simplex/sparse_matrix_ops.h +++ b/src/math/simplex/sparse_matrix_ops.h @@ -26,73 +26,44 @@ namespace simplex { class sparse_matrix_ops { public: static void kernel(sparse_matrix& M, vector>& K) { - mpq_ext::numeral coeff; rational D; vector d, c; - unsigned m = M.num_vars(); + unsigned n = M.num_vars(), m = M.num_rows(); auto& mgr = M.get_manager(); - for (unsigned v = 0; v < m; ++v) - c.push_back(0); + c.resize(m, 0u); + d.resize(n, 0u); - unsigned nullity = 0, rank = 0; - for (auto const& row : M.get_rows()) { - // scan for non-zero variable in row - bool found = false; - d.push_back(0); - ++nullity; - for (auto& [coeff1, v] : M.get_row(row)) { - if (mpq_manager::is_zero(coeff1)) + 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; - if (c[v] != 0) - continue; - --nullity; - ++rank; - d.back() = v + 1; - c[v] = row.id() + 1; - D = rational(-1) / coeff1; - mgr.set(coeff1, mpq(-1)); - // eliminate v from other rows. - for (auto [row2, row_entry2] : M.get_rows(v)) { - if (row.id() >= row2.id()) - continue; - mpq & m_js = row_entry2->m_coeff; - mgr.set(m_js, (D * m_js).to_mpq()); - } - for (auto& [m_ik, w] : M.get_row(row)) { - if (v == w) - continue; - D = m_ik; - mgr.set(m_ik, mpq(0)); - for (auto [row2, row_entry2] : M.get_rows(w)) { - if (row.id() >= row2.id()) - continue; - auto& m_js = M.get_coeff(row2, v); - auto & m_is = row_entry2->m_coeff; - mgr.set(m_is, (m_is + D * m_js).to_mpq()); - } - } + 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; - } + } } - mgr.del(coeff); - - std::cout << "nullity " << nullity << "\n"; - std::cout << "rank " << rank << "\n"; - unsigned_vector pivots(rank, 0u); - unsigned_vector nonpivots(nullity, 0u); - - - - // TODO: extract kernel using d - for (unsigned k = 0; k < d.size(); ++k) { + for (unsigned k = 0; k < n; ++k) { if (d[k] != 0) continue; K.push_back(vector()); - for (unsigned i = 0; i < d.size(); ++i) { + for (unsigned i = 0; i < n; ++i) { if (d[i] > 0) { - // row r = row(i); - // K.back().push_back(M[d[i]-1][k]); + 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)); @@ -102,6 +73,7 @@ namespace simplex { } } + }; diff --git a/src/test/simplex.cpp b/src/test/simplex.cpp index a3be5d3bb..5c52642bd 100644 --- a/src/test/simplex.cpp +++ b/src/test/simplex.cpp @@ -31,23 +31,23 @@ static vector vec(int i, int j, int k) { return nv; } -// static vector vec(int i, int j, int k, int l) { -// vector nv = vec(i, j, k); -// nv.push_back(R(l)); -// return nv; -// } +static vector vec(int i, int j, int k, int l) { + vector nv = vec(i, j, k); + nv.push_back(R(l)); + return nv; +} -/// static vector vec(int i, int j, int k, int l, int x) { -/// vector nv = vec(i, j, k, l); -/// nv.push_back(R(x)); -/// return nv; -/// } +static vector vec(int i, int j, int k, int l, int x) { + vector nv = vec(i, j, k, l); + nv.push_back(R(x)); + return nv; +} -// static vector vec(int i, int j, int k, int l, int x, int y) { -// vector nv = vec(i, j, k, l, x); -// nv.push_back(R(y)); -// return nv; -// } +static vector vec(int i, int j, int k, int l, int x, int y) { + vector nv = vec(i, j, k, l, x); + nv.push_back(R(y)); + return nv; +} // static vector vec(int i, int j, int k, int l, int x, int y, int z) { // vector nv = vec(i, j, k, l, x, y); @@ -133,7 +133,7 @@ static void test4() { } static void add(qmatrix& m, vector const& v) { - m.ensure_var(v.size()); + m.ensure_var(v.size()-1); auto r = m.mk_row(); for (unsigned u = 0; u < v.size(); ++u) m.add_var(r, v[u].to_mpq(), u); @@ -142,12 +142,16 @@ static void add(qmatrix& m, vector const& v) { static void test5() { unsynch_mpq_manager m; qmatrix M(m); - add(M, vec(1, 2, 3)); - add(M, vec(2, 2, 4)); + add(M, vec(1, 0, -3, 0, 2, -8)); + add(M, vec(0, 1, 5, 0, -1, 4)); + add(M, vec(0, 0, 0, 1, 7, -9)); + add(M, vec(0, 0, 0, 0, 0, 0)); M.display(std::cout); vector> K; kernel(M, K); std::cout << "after\n"; + for (auto const& v : K) + std::cout << v << "\n"; M.display(std::cout); }