diff --git a/src/sat/smt/CMakeLists.txt b/src/sat/smt/CMakeLists.txt
index 4ba80ee84..c5bf5b60d 100644
--- a/src/sat/smt/CMakeLists.txt
+++ b/src/sat/smt/CMakeLists.txt
@@ -44,6 +44,7 @@ z3_add_component(sat_smt
     sat_th.cpp
     tseitin_theory_checker.cpp
     user_solver.cpp
+    xor_matrix_finder.cpp
     xor_solver.cpp
   COMPONENT_DEPENDENCIES
     sat
diff --git a/src/sat/smt/xor_matrix_finder.cpp b/src/sat/smt/xor_matrix_finder.cpp
new file mode 100644
index 000000000..f230c7931
--- /dev/null
+++ b/src/sat/smt/xor_matrix_finder.cpp
@@ -0,0 +1,208 @@
+#include "sat/smt/xor_matrix_finder.h"
+#include "sat/smt/xor_solver.h"
+
+
+namespace xr {
+
+    xor_matrix_finder::xor_matrix_finder(solver& s) : m_xor(s), m_sat(s.s()) { }
+        
+    inline bool xor_matrix_finder::belong_same_matrix(const constraint& x) {
+        uint32_t comp_num = -1;
+        for (sat::bool_var v : x) {
+            if (m_table[v] == l_undef) 
+                //Belongs to none, abort
+                return false;            
+            if (comp_num == -1) 
+                //Belongs to this one
+                comp_num = m_table[v];            
+            else if (comp_num != m_table[v]) 
+                //Another var in this XOR belongs to another component
+                return false;            
+        }
+        return true;
+    }
+    
+    bool xor_matrix_finder::find_matrices(bool& can_detach) {
+
+        SASSERT(!m_sat.inconsistent());
+#if 0
+        SASSERT(m_solver->gmatrices.empty());
+        can_detach = true;
+    
+        m_table.clear();
+        m_table.resize(m_solver->nVars(), l_undef);
+        m_reverseTable.clear();
+        clash_vars_unused.clear();
+        m_matrix_no = 0;
+        
+        XorFinder finder(NULL, solver);
+    
+        for(auto& x: m_solver->xorclauses_unused) m_solver->xorclauses.push_back(std::move(x));
+        m_solver->xorclauses_unused.clear();
+        m_solver->clauseCleaner->clean_xor_clauses(solver->xorclauses);
+    
+        finder.grab_mem();
+        finder.move_xors_without_connecting_vars_to_unused();
+        if (!finder.xor_together_xors(m_solver->xorclauses)) return false;
+    
+        finder.move_xors_without_connecting_vars_to_unused();
+        finder.clean_equivalent_xors(m_solver->xorclauses);
+        for(const auto& x: m_solver->xorclauses_unused)
+            clash_vars_unused.insert(x.clash_vars.begin(), x.clash_vars.end());
+    
+        if (m_solver->xorclauses.size() < m_solver->conf.gaussconf.min_gauss_xor_clauses) {
+            can_detach = false;
+            return true;
+        }
+    
+        //Just one giant matrix.
+        if (!m_solver->conf.gaussconf.doMatrixFind) {
+            m_solver->gmatrices.push_back(new EGaussian(m_solver, 0, m_solver->xorclauses));
+            m_solver->gqueuedata.resize(m_solver->gmatrices.size());
+            return true;
+        }
+    
+        std::vector<uint32_t> newSet;
+        set<uint32_t> tomerge;
+        for (const constraint& x : m_solver->m_constraints) {
+            if (belong_same_matrix(x))
+                continue;
+    
+            tomerge.clear();
+            newSet.clear();
+            for (uint32_t v : x) {
+                if (m_table[v] != l_undef)
+                    tomerge.insert(m_table[v]);
+                else
+                    newSet.push_back(v);
+            }
+            if (tomerge.size() == 1) {
+                const uint32_t into = *tomerge.begin();
+                auto intoReverse = m_reverseTable.find(into);
+                for (uint32_t i = 0; i < newSet.size(); i++) {
+                    intoReverse->second.push_back(newSet[i]);
+                    m_table[newSet[i]] = into;
+                }
+                continue;
+            }
+    
+            for (uint32_t v: tomerge) {
+                newSet.insert(newSet.end(), m_reverseTable[v].begin(), m_reverseTable[v].end());
+                m_reverseTable.erase(v);
+            }
+            for (uint32_t i = 0; i < newSet.size(); i++)
+                m_table[newSet[i]] = m_matrix_no;
+            m_reverseTable[m_matrix_no] = newSet;
+            m_matrix_no++;
+        }
+    
+        uint32_t numMatrixes = set_matrixes();
+    
+        const bool time_out =  false;
+        
+        return !m_solver->inconsistent();
+#endif
+        return false;
+    }
+    
+    uint32_t xor_matrix_finder::set_matrixes() {
+        return 0;
+
+        svector<matrix_shape> matrix_shapes;
+        svector<ptr_vector<constraint>> xors_in_matrix(m_matrix_no);
+
+        for (unsigned i = 0; i < m_matrix_no; i++) {
+            matrix_shapes.push_back(matrix_shape(i));
+            matrix_shapes[i].m_num = i;
+            matrix_shapes[i].m_cols = m_reverseTable[i].size();
+        }
+#if 0        
+        for (constraint* x : m_xor.m_constraints) {
+            //take 1st variable to check which matrix it's in.
+            const uint32_t matrix = m_table[(*x)[0]];
+            SASSERT(matrix < m_matrix_no);
+    
+            //for stats
+            matrix_shapes[matrix].m_rows ++;
+            matrix_shapes[matrix].m_sum_xor_sizes += x->get_size();
+            xors_in_matrix[matrix].push_back(x);
+        }
+        m_solver->m_constraints.clear();
+    
+        for(auto& m: matrix_shapes) {
+            if (m.tot_size() > 0) {
+                m.m_density = (double)m.m_sum_xor_sizes / (double)(m.tot_size());
+            }
+        }
+    
+        std::sort(matrix_shapes.begin(), matrix_shapes.end(), m_sorter);
+    
+        uint32_t realMatrixNum = 0;
+        uint32_t unusedMatrix = 0;
+        uint32_t too_few_rows_matrix = 0;
+        uint32_t unused_matrix_printed = 0;
+        for (int a = m_matrix_no - 1; a >= 0; a--) {
+            matrix_shape& m = matrix_shapes[a];
+            uint32_t i = m.m_num;
+            if (m.m_rows == 0) {
+                continue;
+            }
+    
+            bool use_matrix = true;
+    
+    
+            //Over- or undersized
+            if (use_matrix && m.m_rows > m_solver->conf.gaussconf.max_matrix_rows) {
+                use_matrix = false;
+                // Too many rows in matrix
+            }
+            if (use_matrix && m.m_cols > m_solver->conf.gaussconf.max_matrix_columns) {
+                use_matrix = false;
+                // Too many columns in matrix
+            }
+    
+            if (use_matrix && m.m_rows < m_solver->conf.gaussconf.min_matrix_rows) {
+                use_matrix = false;
+                too_few_rows_matrix++;
+                // Too few rows in matrix
+            }
+    
+            //Over the max number of matrixes
+            if (use_matrix && realMatrixNum >= m_solver->conf.gaussconf.max_num_matrices) {
+                // above max number of matrixes
+                use_matrix = false;
+            }
+    
+            //if already detached, we MUST use the matrix
+            for(const auto& x: xors_in_matrix[i]) {
+                if (x->is_detached()) {
+                    use_matrix = true;
+                    break;
+                }
+            }
+    
+            if (m_solver->conf.force_use_all_matrixes) {
+                use_matrix = true;
+            }
+    
+            if (use_matrix) {
+                m_solver->gmatrices.push_back(
+                    new EGaussian(m_solver, realMatrixNum, xors_in_matrix[i]));
+                m_solver->gqueuedata.resize(m_solver->gmatrices.size());
+    
+                realMatrixNum++;
+                SASSERT(m_solver->gmatrices.size() == realMatrixNum);
+            } else {
+                for (auto& x: xors_in_matrix[i]) {
+                    m_solver->xorclauses_unused.push_back(x);
+                    //cout<< "c [matrix]xor not in matrix, now unused_xors size: " << unused_xors.size() << endl;
+                    clash_vars_unused.insert(x.clash_vars.begin(), x.clash_vars.end());
+                }
+                unusedMatrix++;
+            }
+        }
+    
+        return realMatrixNum;
+#endif
+    }
+}
diff --git a/src/sat/smt/xor_matrix_finder.h b/src/sat/smt/xor_matrix_finder.h
new file mode 100644
index 000000000..bbec3424c
--- /dev/null
+++ b/src/sat/smt/xor_matrix_finder.h
@@ -0,0 +1,57 @@
+#pragma once
+
+#include "util/debug.h"
+#include "util/vector.h"
+#include "util/uint_set.h"
+#include "util/map.h"
+#include "sat/sat_solver.h"
+
+namespace xr {
+
+    class solver;
+    class constraint;
+    
+    class xor_matrix_finder {
+
+                        
+        struct matrix_shape {
+            matrix_shape(uint32_t matrix_num) : m_num(matrix_num) {}
+        
+            matrix_shape() {}
+        
+            uint32_t m_num = 0;
+            uint32_t m_rows = 0;
+            uint32_t m_cols = 0;
+            uint32_t m_sum_xor_sizes = 0;
+            double m_density = 0;
+        
+            uint64_t tot_size() const {
+                return (uint64_t)m_rows*(uint64_t)m_cols;
+            }
+        };
+
+        struct sorter {
+            bool operator () (const matrix_shape& left, const matrix_shape& right) {
+                return left.m_sum_xor_sizes < right.m_sum_xor_sizes;
+            }
+        };
+        
+        u_map<unsigned_vector> m_reverseTable; //matrix -> vars
+        unsigned_vector        m_table; //var -> matrix
+        unsigned               m_matrix_no = 0;
+        sorter                 m_sorter;
+        solver&                m_xor;
+        sat::solver&           m_sat;
+        
+        
+        unsigned set_matrixes();
+        
+        inline bool belong_same_matrix(const constraint& x);
+        
+    public:
+        xor_matrix_finder(solver& solver);
+        
+        bool find_matrices(bool& can_detach);
+        uint_set clash_vars_unused;
+    };
+}
diff --git a/src/sat/smt/xor_solver.h b/src/sat/smt/xor_solver.h
index eb88e0717..6613c3a6d 100644
--- a/src/sat/smt/xor_solver.h
+++ b/src/sat/smt/xor_solver.h
@@ -17,6 +17,38 @@ Abstract:
 #include "sat/smt/euf_solver.h"
 
 namespace xr {
+
+    class constraint {
+        size_t          m_size;
+        bool            m_detached;
+        size_t          m_obj_size;
+        bool            m_rhs;
+        sat::bool_var   m_vars[0];
+        
+    public:
+        static size_t get_obj_size(unsigned num_lits) { return sat::constraint_base::obj_size(sizeof(constraint) + num_lits * sizeof(sat::bool_var)); }
+        
+        constraint(const svector<sat::bool_var>& ids, bool expected_result) : m_size(ids.size()), m_detached(false), m_obj_size(get_obj_size(ids.size())), m_rhs(expected_result) {
+            unsigned i = 0;
+            for (auto v : ids)
+                m_vars[i++] = v;
+        }
+        sat::ext_constraint_idx cindex() const { return sat::constraint_base::mem2base(this); }
+        void deallocate(small_object_allocator& a) { a.deallocate(m_obj_size, sat::constraint_base::mem2base_ptr(this)); }
+        sat::bool_var operator[](unsigned i) const { return m_vars[i]; }
+        bool is_detached() const { return m_detached; }
+        size_t get_size() const { return m_size; }
+        bool get_rhs() const { return m_rhs; }
+        sat::bool_var const* begin() const { return m_vars; }
+        sat::bool_var const* end() const { return m_vars + m_size; }
+        std::ostream& display(std::ostream& out) const {
+            bool first = true;
+            for (auto v : *this)
+                out << (first ? "" : " ^ ") << v, first = false;
+            return out << " = " << m_rhs;
+        }
+    };
+    
     class solver : public euf::th_solver {
         euf::solver* m_ctx = nullptr;
         sat::sat_internalizer& si;