Symbolic gaussian elimination prototype (#128)

Components of this prototype have been incorporated in the kinetic schemes PR #145, but that PR does not retain the more extensive numeric validation code nor does it perform the operation or fill in accounting that is reported in the demo for the evaluation of various pivot choice schemes.

I propose that we keep this demo application for its utility and code base when we revisit this problem, particularly for the investigation of:

the numerical behaviour for larger matrices (limited dynamic range will likely force periodic normalization of temporaries);
better pivot selection schemes for reduced fill-in;
symbolic computation that will arise in higher order implicit schemes.

CMakeLists.txt

@@ -186,6 +186,7 @@ endif()
 
 add_subdirectory(mechanisms)
 add_subdirectory(src)
+add_subdirectory(proto/symge)
 add_subdirectory(tests)
 add_subdirectory(miniapp)
 

proto/symge/CMakeLists.txt

+add_executable(symge-demo symge-demo.cpp)
+
+set_target_properties(symge-demo
+   PROPERTIES
+   RUNTIME_OUTPUT_DIRECTORY "${CMAKE_BINARY_DIR}/proto"
+)
+

proto/symge/msparse.hpp

+#pragma once
+
+#include <algorithm>
+#include <utility>
+#include <initializer_list>
+#include <iterator>
+#include <vector>
+
+#include <util/compat.hpp>
+#include <util/iterutil.hpp>
+#include <util/rangeutil.hpp>
+
+namespace msparse {
+
+namespace util = nest::mc::util;
+
+struct msparse_error: std::runtime_error {
+    msparse_error(const std::string &what): std::runtime_error(what) {}
+};
+
+constexpr unsigned row_npos = unsigned(-1);
+
+template <typename X>
+class row {
+public:
+    struct entry {
+        unsigned col;
+        X value;
+    };
+    static constexpr unsigned npos = row_npos;
+
+private:
+    // entries must have strictly monotonically increasing col numbers.
+    std::vector<entry> data_;
+
+    bool check_invariant() const {
+        for (unsigned i = 1; i<data_.size(); ++i) {
+            if (data_[i].col<=data_[i-1].col) return false;
+        }
+        return true;
+    }
+
+public:
+    row() = default;
+    row(const row&) = default;
+
+    row(std::initializer_list<entry> il): data_(il) {
+        if (!check_invariant())
+            throw msparse_error("improper row element list");
+    }
+
+    template <typename InIter>
+    row(InIter b, InIter e): data_(b, e) {
+        if (!check_invariant())
+            throw msparse_error("improper row element list");
+    }
+
+    unsigned size() const { return data_.size(); }
+    bool empty() const { return size()==0; }
+
+    auto begin() -> decltype(data_.begin()) { return data_.begin(); }
+    auto begin() const -> decltype(data_.cbegin()) { return data_.cbegin(); }
+    auto end() -> decltype(data_.end()) { return data_.end(); }
+    auto end() const -> decltype(data_.cend()) { return data_.cend(); }
+
+    unsigned mincol() const {
+        return empty()? npos: data_.front().col;
+    }
+
+    unsigned mincol_after(unsigned c) const {
+        auto i = std::upper_bound(data_.begin(), data_.end(), c,
+            [](unsigned a, const entry& b) { return a<b.col; });
+
+        return i==data_.end()? npos: i->col;
+    }
+
+    unsigned maxcol() const {
+        return empty()? npos: data_.back().col;
+    }
+
+    const entry& get(unsigned i) const {
+        return data_[i];
+    }
+
+    void push_back(const entry& e) {
+        if (!empty() && e.col <= data_.back().col)
+            throw msparse_error("cannot push_back row elements out of order");
+        data_.push_back(e);
+    }
+
+    unsigned index(unsigned c) const {
+        auto i = std::lower_bound(data_.begin(), data_.end(), c,
+            [](const entry& a, unsigned b) { return a.col<b; });
+
+        return (i==data_.end() || i->col!=c)? npos: std::distance(data_.begin(), i);
+    }
+
+    // remove all entries from column c onwards
+    void truncate(unsigned c) {
+        auto i = std::lower_bound(data_.begin(), data_.end(), c,
+            [](const entry& a, unsigned b) { return a.col<b; });
+        data_.erase(i, data_.end());
+    }
+
+    X operator[](unsigned c) const {
+        auto i = index(c);
+        return i==npos? X{}: data_[i].value;
+    }
+
+    struct assign_proxy {
+        row<X>& row_;
+        unsigned c;
+
+        assign_proxy(row<X>& r, unsigned c): row_(r), c(c) {}
+
+        operator X() const { return const_cast<const row<X>&>(row_)[c]; }
+        assign_proxy& operator=(const X& x) {
+            auto i = std::lower_bound(row_.data_.begin(), row_.data_.end(), c,
+                [](const entry& a, unsigned b) { return a.col<b; });
+
+            if (i==row_.data_.end() || i->col!=c) {
+                row_.data_.insert(i, {c, x});
+            }
+            else if (x == X{}) {
+                row_.data_.erase(i);
+            }
+            else {
+                i->value = x;
+            }
+
+            return *this;
+        }
+    };
+
+    assign_proxy operator[](unsigned c) {
+        return assign_proxy{*this, c};
+    }
+
+    template <typename RASeq>
+    auto dot(const RASeq& v) const -> decltype(X{}*util::front(v)) {
+        using result_type = decltype(X{}*util::front(v));
+        result_type s{};
+
+        auto nv = util::size(v);
+        for (const auto& e: data_) {
+            if (e.col>=nv) throw msparse_error("right multiplicand too short");
+            s += e.value*v[e.col];
+        }
+        return s;
+    }
+};
+
+template <typename X>
+class matrix {
+private:
+    std::vector<row<X>> rows;
+    unsigned cols = 0;
+    unsigned aug = row_npos;
+
+public:
+    static constexpr unsigned npos = row_npos;
+
+    matrix() = default;
+    matrix(unsigned n, unsigned c): rows(n), cols(c) {}
+
+    row<X>& operator[](unsigned i) { return rows[i]; }
+    const row<X>& operator[](unsigned i) const { return rows[i]; }
+
+    unsigned size() const { return rows.size(); }
+    unsigned nrow() const { return size(); }
+    unsigned ncol() const { return cols; }
+    unsigned augcol() const { return aug; }
+
+    bool empty() const { return size()==0; }
+    bool augmented() const { return aug!=npos; }
+
+    template <typename Seq>
+    void augment(const Seq& col_dense) {
+        unsigned r = 0;
+        for (const auto& v: col_dense) {
+            if (r>=rows.size()) throw msparse_error("augmented column size mismatch");
+            rows[r++].push_back({cols, v});
+        }
+        if (aug==npos) aug=cols;
+        ++cols;
+    }
+
+    void diminish() {
+        if (aug==npos) return;
+        for (auto& row: rows) row.truncate(aug);
+        cols = aug;
+        aug = npos;
+    }
+};
+
+// sparse * dense vector muliply: writes A*x to b.
+template <typename AT, typename RASeqX, typename SeqB>
+void mul_dense(const matrix<AT>& A, const RASeqX& x, SeqB& b) {
+    auto bi = std::begin(b);
+    unsigned n = A.nrow();
+    for (unsigned i = 0; i<n; ++i) {
+        if (bi==compat::end(b)) throw msparse_error("output sequence b too short");
+        *bi++ = A[i].dot(x);
+    }
+}
+
+} // namespace msparse

proto/symge/symbolic.hpp

+#pragma once
+
+#include <iosfwd>
+#include <string>
+#include <stdexcept>
+#include <vector>
+
+#include <util/optional.hpp>
+
+namespace sym {
+
+template <typename X>
+using optional = nest::mc::util::optional<X>;
+
+using nest::mc::util::nothing;
+using nest::mc::util::just;
+
+struct symbol_error: public std::runtime_error {
+    symbol_error(const std::string& what): std::runtime_error(what) {}
+};
+
+// Note impl definitions below over template argument S are used to resolve
+// class declaration dependencies.
+
+namespace impl {
+    // Represents product of two symbols, or zero if either symbol is null.
+    template <typename S>
+    struct symbol_term_ {
+        S a, b;
+
+        symbol_term_() = default;
+
+        // true iff representing non-zero
+        operator bool() const {
+            return a && b;
+        }
+    };
+
+    // Represents the difference between two product terms.
+    template <typename S>
+    struct symbol_term_diff_ {
+        symbol_term_<S> left, right;
+
+        symbol_term_diff_() = default;
+        symbol_term_diff_(const symbol_term_<S>& left): left(left), right{} {}
+        symbol_term_diff_(const symbol_term_<S>& left, const symbol_term_<S>& right):
+           left(left), right(right) {}
+    };
+
+    // A symbol can be primitive (has no expansion) or is defined as a product term difference.
+    template <typename S>
+    using symbol_def_ = optional<symbol_term_diff_<S>>;
+
+    // A symbol table represents a set of symbol names and definitions, referenced
+    // by symbol index.
+    template <typename S>
+    class symbol_table_ {
+    public:
+        struct table_entry {
+            std::string name;
+            symbol_def_<S> def;
+        };
+
+        S define(const std::string& name, const symbol_def_<S>& definition = nothing) {
+            unsigned idx = size();
+            entries.push_back({name, definition});
+            return S{idx, this};
+        }
+
+        S operator[](unsigned i) const {
+            if (i>=size()) throw symbol_error("no such symbol");
+            return S{i, this};
+        }
+
+        std::size_t size() const {
+            return entries.size();
+        }
+
+        const symbol_def_<S>& def(S s) const {
+            if (!valid(s)) throw symbol_error("symbol not present in this table");
+            return entries[s.idx].def;
+        }
+
+        const std::string& name(S s) const {
+            if (!valid(s)) throw symbol_error("symbol not present in this table");
+            return entries[s.idx].name;
+        }
+
+    private:
+        std::vector<table_entry> entries;
+        bool valid(S s) const {
+            return s.tbl==this && s.idx<entries.size();
+        }
+    };
+} // namespace impl
+
+class symbol {
+public:
+    using symbol_table = impl::symbol_table_<symbol>;
+    using symbol_def = impl::symbol_def_<symbol>;
+
+private:
+    friend class impl::symbol_table_<symbol>;
+    unsigned idx = 0;
+    const symbol_table* tbl = nullptr;
+
+    // Symbols are created through a symbol table, or are
+    // default-constructed 'null' symbols.
+    symbol(unsigned idx, const symbol_table* tbl):
+        idx(idx), tbl(tbl) {}
+
+public:
+    symbol() = default;
+    symbol(const symbol&) = default;
+    symbol& operator=(const symbol&) = default;
+
+    // A symbol is 'null' if it has no corresponding symbol table.
+    operator bool() const { return (bool)tbl; }
+
+    std::string str() const {
+        return tbl? tbl->name(*this): "";
+    }
+
+    symbol_def def() const {
+        return tbl? tbl->def(*this): symbol_def{};
+    }
+
+    bool primitive() const { return (bool)def(); }
+
+    optional<unsigned> index(const symbol_table& in_table) const {
+        return tbl==&in_table? just(idx): nothing;
+    }
+};
+
+using symbol_term = impl::symbol_term_<symbol>;
+using symbol_term_diff = impl::symbol_term_diff_<symbol>;
+using symbol_def = impl::symbol_def_<symbol>;
+using symbol_table = impl::symbol_table_<symbol>;
+
+inline symbol_term_diff operator-(const symbol_term& left, const symbol_term& right) {
+    return symbol_term_diff{left, right};
+}
+
+inline symbol_term_diff operator-(const symbol_term& right) {
+    return symbol_term_diff{symbol_term{}, right};
+}
+
+inline symbol_term operator*(const symbol& a, const symbol& b) {
+    return symbol_term{a, b};
+}
+
+inline std::ostream& operator<<(std::ostream& o, const symbol& s) {
+    return o << s.str();
+}
+
+inline std::ostream& operator<<(std::ostream& o, const symbol_term& term) {
+    if (term) return o << term.a.str() << '*' << term.b.str();
+    else return o << '0';
+}
+
+inline std::ostream& operator<<(std::ostream& o, const symbol_term_diff& diff) {
+    if (!diff.right) return o << diff.left;
+    else {
+        if (diff.left) o << diff.left;
+        o << '-';
+        return o << diff.right;
+    }
+}
+
+// A store represents a map from symbols (all from one table) to values.
+class store {
+private:
+    const symbol_table& table;
+    std::vector<optional<double>> data;
+
+public:
+    explicit store(const symbol_table& table): table(table) {}
+
+    optional<double>& operator[](const symbol& s) {
+        if (auto idx = s.index(table)) {
+            if (*idx>=data.size()) data.resize(1+*idx);
+            return data[*idx];
+        }
+        throw symbol_error("symbol not associated with store table");
+    }
+
+    optional<double> operator[](const symbol& s) const {
+        auto idx = s.index(table);
+        if (idx && *idx<data.size()) return data[*idx];
+        else return nothing;
+    }
+
+    optional<double> evaluate(const symbol& s) {
+        auto& value = (*this)[s];
+        if (!value) {
+            if (auto def = s.def()) {
+                value = evaluate(*def);
+            }
+        }
+        return value;
+    }
+
+    optional<double> evaluate(const symbol_term& t) {
+        if (!t) return 0;
+        auto a = evaluate(t.a);
+        auto b = evaluate(t.b);
+        return a && b? ++mul_count, just((*a)*(*b)): nothing;
+    }
+
+    optional<double> evaluate(const symbol_term_diff& d) {
+        auto l = evaluate(d.left);
+        auto r = evaluate(d.right);
+        return l && r? ++sub_count, just((*l)-(*r)): nothing;
+    }
+
+    // stat collection...
+    unsigned mul_count = 0;
+    unsigned sub_count = 0;
+};
+
+} // namespace sym
+

proto/symge/symge-demo.cpp

+// bloody CMake
+#ifdef NDEBUG
+#undef NDEBUG
+#endif
+
+#include <cassert>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+#include <queue>
+#include <random>
+#include <sstream>
+#include <string>
+#include <unordered_set>
+#include <utility>
+
+#include <util/optional.hpp>
+
+#include "symbolic.hpp"
+#include "msparse.hpp"
+
+using namespace sym;
+
+// Identifier name picking helper routines
+
+void join_impl(std::ostream& ss, const std::string& sep) {}
+
+template <typename Head, typename... Args>
+void join_impl(std::ostream& ss, const std::string& sep, Head&& head, Args&&... tail) {
+    if (sizeof...(tail)==0)
+        ss << std::forward<Head>(head);
+    else
+        join_impl(ss << std::forward<Head>(head) << sep, sep, std::forward<Args>(tail)...);
+}
+
+// Return string representation of all arguments, separated by `sep`.
+template <typename... Args>
+std::string join(const std::string& sep, Args&&... items) {
+    std::stringstream ss;
+    join_impl(ss, sep, std::forward<Args>(items)...);
+    return ss.str();
+}
+
+// Maintain a collection of unique identifiers.
+class id_maker {
+private:
+    std::unordered_set<std::string> ids;
+
+public:
+    // Find the next string lexicographically after argument, allowing only
+    // ASCII letters and digits to change.
+    static std::string next_id(std::string s) {
+        unsigned l = s.size();
+        if (l==0) return "a";
+
+        unsigned i = l-1;
+
+        constexpr char l0='a', l1='z';
+        constexpr char u0='A', u1='Z';
+        constexpr char d0='0', d1='9';
+        for (;;) {
+            char& c = s[i];
+            if ((c>=l0 && c<l1) || (c>=u0 && c<u1) || (c>=d0 && c<d1)) {
+                ++c;
+                return s;
+            }
+            if (c==l1) c=l0;
+            if (c==u1) c=u0;
+            if (c==d1) c=d0;
+            if (i==0) break;
+            --i;
+        }
+
+        // prepend a character based on class of first
+        if (s[0]==u0) return u0+s;
+        if (s[0]==d0) return d0+s;
+        return l0+s;
+    }
+
+    // Make a new identifier based on the given arguments.
+    template <typename... Args>
+    std::string operator()(Args&&... elements) {
+        std::string name = join("",std::forward<Args>(elements)...);
+        if (name.empty()) name = "a";
+
+        while (ids.count(name)) {
+            name = next_id(name);
+        }
+        ids.insert(name);
+        return name;
+    }
+
+    void reserve(std::string name) {
+        ids.insert(std::move(name));
+    }
+};
+
+// Output helper functions
+
+template <typename X>
+std::ostream& operator<<(std::ostream& o, const optional<X>& x) {
+    return x? o << *x: o << "nothing";
+}
+
+template <typename X>
+std::ostream& operator<<(std::ostream& o, const msparse::matrix<X>& m) {
+    for (unsigned r = 0; r<m.nrow(); ++r) {
+        o << '|';
+        for (unsigned c = 0; c<m.ncol(); ++c) {
+            if (c==m.augcol()) o << " | ";
+            o << std::setw(12) << m[r][c];
+        }
+        o << " |\n";
+    }
+    return o;
+}
+
+template <typename Sep, typename V>
+struct sepval_proxy {
+    const Sep& sep;
+    const V& v;
+
+    sepval_proxy(const Sep& sep, const V& v): sep(sep), v(v) {}
+
+    friend std::ostream& operator<<(std::ostream& O, const sepval_proxy& sv) {
+        bool first = true;
+        for (const auto& x: sv.v) {
+            if (!first) O << sv.sep;
+            first = false;
+            O << x;
+        }
+        return O;
+    }
+};
+
+template <typename Sep, typename V>
+sepval_proxy<Sep,V> sepval(const Sep& sep, const V& v) { return sepval_proxy<Sep,V>(sep, v); }
+
+std::ostream& operator<<(std::ostream& o, const symbol_table& syms) {
+    for (unsigned i = 0; i<syms.size(); ++i) {
+        symbol s = syms[i];
+        if (s.def()) o << s << ": " << s.def() << "\n";
+    }
+    return o;
+}
+
+// Symbolic GE
+
+using sym_row = msparse::row<symbol>;
+using sym_matrix = msparse::matrix<symbol>;
+
+// Returns q[c]*p - p[c]*q
+template <typename DefineSym>
+sym_row row_reduce(unsigned c, const sym_row& p, const sym_row& q, DefineSym& define_sym) {
+    if (p.index(c)==p.npos || q.index(c)==q.npos) throw std::runtime_error("improper row GE");
+
+    sym_row u;
+    symbol x = q[c];
+    symbol y = p[c];
+
+    auto piter = p.begin();
+    auto qiter = q.begin();
+    unsigned pj = piter->col;
+    unsigned qj = qiter->col;
+
+    while (piter!=p.end() || qiter!=q.end()) {
+        unsigned j = std::min(pj, qj);
+        symbol_term t1, t2;
+
+        if (j==pj) {
+            t1 = x*piter->value;
+            ++piter;
+            pj = piter==p.end()? p.npos: piter->col;
+        }
+        if (j==qj) {
+            t2 = y*qiter->value;
+            ++qiter;
+            qj = qiter==q.end()? q.npos: qiter->col;
+        }
+        if (j!=c) {
+            u.push_back({j, define_sym(t1-t2)});
+        }
+    }
+    return u;
+}
+
+// Actual GE reduction.
+template <typename DefineSym>
+void gj_reduce_simple(sym_matrix& A, unsigned ncol, const DefineSym& define_sym) {
+    // Expect A to be a (possibly column-augmented) diagonally dominant
+    // matrix; take diagonal elements as pivots.
+    if (A.nrow()>A.ncol()) throw std::runtime_error("improper matrix for reduction");
+
+    for (unsigned pivrow = 0; pivrow<A.nrow(); ++pivrow) {
+        unsigned pivcol = pivrow;
+
+        for (unsigned i = 0; i<A.nrow(); ++i) {
+            if (i==pivrow || A[i].index(pivcol)==msparse::row_npos) continue;
+
+            A[i] = row_reduce(pivcol, A[i], A[pivrow], define_sym);
+        }
+    }
+}
+
+// simple greedy estimate based on immediate fill cost
+double estimate_cost(const sym_matrix& A, unsigned p) {
+    struct count_sym_mul {
+        mutable unsigned nmul = 0;
+        mutable unsigned nfill = 0;
+        symbol operator()(symbol_term_diff t) {
+            bool l = t.left, r = t.right;
+            nmul += l+r;
+            nfill += r&!l;
+            return symbol{};
+        }
+    };
+
+    count_sym_mul counter;
+    for (unsigned i = 0; i<A.nrow(); ++i) {
+        if (i==p || A[i].index(p)==msparse::row_npos) continue;
+
+        row_reduce(p, A[i], A[p], counter);
+    }
+    return counter.nfill;
+}
+
+template <typename DefineSym>
+void gj_reduce(sym_matrix& A, unsigned ncol, const DefineSym define_sym) {
+    // Expect A to be a (possibly column-augmented) diagonally dominant
+    // matrix; take diagonal elements as pivots.
+    if (A.nrow()>A.ncol()) throw std::runtime_error("improper matrix for reduction");
+
+    std::vector<unsigned> pivots;
+    for (unsigned r = 0; r<A.nrow(); ++r) {
+        pivots.push_back(r);
+    }
+
+    std::vector<double> cost(pivots.size());
+
+    while (!pivots.empty()) {
+        for (unsigned i = 0; i<pivots.size(); ++i) {
+            cost[pivots[i]] = estimate_cost(A, pivots[i]);
+        }
+
+        std::sort(pivots.begin(), pivots.end(),
+            [&](unsigned r1, unsigned r2) { return cost[r1]>cost[r2]; });
+
+        unsigned pivrow = pivots.back();
+        pivots.erase(std::prev(pivots.end()));
+
+        unsigned pivcol = pivrow;
+
+        for (unsigned i = 0; i<A.nrow(); ++i) {
+            if (i==pivrow || A[i].index(pivcol)==msparse::row_npos) continue;
+
+            A[i] = row_reduce(pivcol, A[i], A[pivrow], define_sym);
+        }
+    }
+}
+
+// Validation
+
+template <typename Rng>
+msparse::matrix<double> make_random_matrix(unsigned n, double row_nnz, Rng& R) {
+    std::uniform_real_distribution<double> U;
+    msparse::matrix<double> M(n, n);
+
+    row_nnz = std::min(row_nnz, n-1.5); // always allow some chance of zeroes.
+    double density = n==1? 1: row_nnz/(n-1);
+    for (unsigned i = 0; i<n; ++i) {
+        for (unsigned j = 0; j<n; ++j) {
+            if (i!=j && U(R)>density) continue;
+            double u = U(R);
+            M[i][j] = i==j? 1+0.2*u: (u-0.5)/n;
+        }
+    }
+    return M;
+}
+
+struct symge_stats {
+    unsigned n;
+    unsigned nnz;
+    unsigned nmul;
+    unsigned nsub;
+    unsigned nsym;
+    double relerr;
+};
+
+template <typename Rng>
+symge_stats run_symge_validation(Rng& R, unsigned n, double row_nnz, bool use_estimator, bool debug = false) {
+    std::uniform_real_distribution<double> U;
+    symge_stats stats = { n, 0, 0, 0, 0, 0. };
+
+    msparse::matrix<double> A = make_random_matrix(n, row_nnz, R);
+    sym_matrix S(n, n);
+
+    symbol_table syms;
+    store values(syms);
+    id_maker make_id;
+
+    // make symbolic matrix with same sparsity pattern as M
+    for (unsigned i = 0; i<A.nrow(); ++i) {
+        sym_row srow;
+        for (const auto& a: A[i]) {
+            unsigned j = a.col;
+            auto s = syms.define(make_id("a", i, j));
+            values[s] = a.value;
+            srow.push_back({j, s});
+            ++stats.nnz;
+        }
+        S[i] = srow;
+    }
+
+    // pick x and compute b = Ax
+    std::vector<double> x(n);
+    for (auto& elem: x) elem = 10*U(R);
+
+    std::vector<double> b(n);
+    mul_dense(A, x, b);
+
+    // augment M and S with rhs
+    A.augment(b);
+    std::vector<symbol> rhs;
+    for (unsigned i = 0; i<n; ++i) {
+        auto s = syms.define(make_id("b", i));
+        rhs.push_back(s);
+        values[s] = b[i];
+    }
+    S.augment(rhs);
+
+    if (debug) {
+        std::cerr << "A|b:\n" << A << "\n";
+        std::cerr << "S:\n" << S << "\n";
+    }
+
+    // perform GE
+    auto nprim_sym = syms.size();
+    auto define_sym = [&](const symbol_def& def) { return syms.define(make_id("t00"), def); };
+    if (use_estimator) {
+        gj_reduce(S, n, define_sym);
+    }
+    else {
+        gj_reduce_simple(S, n, define_sym);
+    }
+    stats.nsym = syms.size()-nprim_sym;
+    if (debug) {
+        std::cerr << "reduced S:\n" << S << "\n";
+        std::cerr << "definitions:\n" << syms << "\n";
+    }
+
+    // validate: compute solution y from reduced S and symbol defs
+    std::vector<double> y(n);
+    double maxerr = 0;
+    double maxx = 0;
+    for (unsigned i = 0; i<n; ++i) {
+        const sym_row& row = S[i];
+        if (row.size()!=2 || row.maxcol()!=n)
+            throw std::runtime_error("unexpected matrix layout!");
+
+        unsigned idx = row.get(0).col;
+        symbol coeff = row.get(0).value;
+        symbol rhs = row.get(1).value;
+
+        if (debug) {
+            std::cerr << "y" << idx << " = " << rhs << "/" << coeff << "\n";
+            std::cerr << "   = " << values.evaluate(rhs).get() << "/" << values.evaluate(coeff).get() << "\n";
+        }
+        y[idx] = values.evaluate(rhs).get()/values.evaluate(coeff).get();
+        maxerr = std::max(maxerr, std::abs(y[idx]-x[idx]));
+        maxx = std::max(maxx, std::abs(x[idx]));
+    }
+    stats.relerr = maxerr/maxx;
+    stats.nsub = values.sub_count;
+    stats.nmul = values.mul_count;
+    if (debug) {
+        std::cerr << "x:\n" << sepval(", ", x) << "\n";
+        std::cerr << "computed solution:\n" << sepval(", ", y) << "\n";
+        std::cerr << "operation count: " << values.mul_count << " mul; " << values.sub_count << " sub\n";
+    }
+    return stats;
+}
+
+void usage() {
+    std::cout << "Usage: symge-demo [-v] [-d] [-e]\n"
+              << "Options:\n"
+              << "  -v      Verbose output: display statistics as table\n"
+              << "  -d      Debug output to stderr\n"
+              << "  -e      Use a simple cost estimator to select pivots\n"
+              << "  -n N    Specify max matrix size (default 10)\n"
+              << "  -k N    Number of matrices to test of each size (default 100)\n"
+              << "  -r R    Mean number of off-diagonal nonzero elements per row (default 2.0)\n";
+}
+
+void usage_error(const std::string& msg) {
+    std::cerr << "symge-demo: " << msg << "\n"
+              << "Try 'symge-demo -h' for more information.\n";
+}
+
+int main(int argc, const char** argv) {
+    bool debug = false;
+    bool verbose = false;
+    bool use_estimator = false;
+    int max_n = 10;
+    int count = 100;
+    double row_mean_nnz = 2.0;
+
+    for (const char** arg=argv+1; arg!=argv+argc; ++arg) {
+        if (!std::strcmp("-v", *arg)) {
+            verbose = true;
+        }
+        else if (!std::strcmp("-d", *arg)) {
+            debug = true;
+        }
+        else if (!std::strcmp("-e", *arg)) {
+            use_estimator = true;
+        }
+        else if (!std::strcmp("-n", *arg)) {
+            if (!arg[1]) {
+                usage_error("expected argument for '-n'");
+                return 2;
+            }
+            max_n = std::atoi(*++arg);
+        }
+        else if (!std::strcmp("-r", *arg)) {
+            if (!arg[1]) {
+                usage_error("expected argument for '-r'");
+                return 2;
+            }
+            row_mean_nnz = std::stod(*++arg);
+        }
+        else if (!std::strcmp("-k", *arg)) {
+            if (!arg[1]) {
+                usage_error("expected argument for '-k'");
+                return 2;
+            }
+            count = std::atoi(*++arg);
+        }
+        else if (!std::strcmp("-h", *arg)) {
+            usage();
+            return 0;
+        }
+        else {
+            std::cerr << "symge-demo: unrecognized argument\n"
+                      << "Try 'symge-demo -h' for more information.\n";
+            return 2;
+        }
+    }
+
+    if (verbose) {
+        char line[80];
+        std::snprintf(line, sizeof(line), "%10s%10s%10s%10s%10s%10s\n",
+            "n", "nnz", "nmul", "nsub", "nsym", "relerr");
+        std::cout << line;
+    }
+
+    std::minstd_rand R;
+    for (unsigned n = 1; n<=max_n; ++n) {
+        for (unsigned k = 1; k<=count; ++k) {
+            auto stats = run_symge_validation(R, n, row_mean_nnz, use_estimator, debug);
+            if (verbose) {
+                char line[80];
+                std::snprintf(line, sizeof(line), "%10d%10d%10d%10d%10d%10lf\n",
+                    stats.n, stats.nnz, stats.nmul, stats.nsub, stats.nsym, stats.relerr);
+
+                std::cout << line;
+            }
+            assert(stats.relerr<1e-6);
+        }
+    }
+
+    return 0;
+}