diff --git a/src/math/lp/CMakeLists.txt b/src/math/lp/CMakeLists.txt
index 97073736717..5d21bba70ba 100644
--- a/src/math/lp/CMakeLists.txt
+++ b/src/math/lp/CMakeLists.txt
@@ -42,7 +42,6 @@ z3_add_component(lp
     random_updater.cpp      
     row_eta_matrix.cpp
     scaler.cpp
-    square_dense_submatrix.cpp
     square_sparse_matrix.cpp
     static_matrix.cpp
   COMPONENT_DEPENDENCIES
diff --git a/src/math/lp/square_dense_submatrix.cpp b/src/math/lp/square_dense_submatrix.cpp
deleted file mode 100644
index 4d9fcec131e..00000000000
--- a/src/math/lp/square_dense_submatrix.cpp
+++ /dev/null
@@ -1,48 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#include <memory>
-#include "util/vector.h"
-#include "math/lp/square_dense_submatrix_def.h"
-template void lp::square_dense_submatrix<double, double>::init(lp::square_sparse_matrix<double, double>*, unsigned int);
-template lp::square_dense_submatrix<double, double>::square_dense_submatrix(lp::square_sparse_matrix<double, double>*, unsigned int);
-template void lp::square_dense_submatrix<double, double>::update_parent_matrix(lp::lp_settings&);
-template bool lp::square_dense_submatrix<double, double>::is_L_matrix() const;
-template void lp::square_dense_submatrix<double, double>::conjugate_by_permutation(lp::permutation_matrix<double, double>&);
-template int lp::square_dense_submatrix<double, double>::find_pivot_column_in_row(unsigned int) const;
-template void lp::square_dense_submatrix<double, double>::pivot(unsigned int, lp::lp_settings&);
-template lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::square_dense_submatrix(lp::square_sparse_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >*, unsigned int);
-template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::update_parent_matrix(lp::lp_settings&);
-template bool lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::is_L_matrix() const;
-template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::conjugate_by_permutation(lp::permutation_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&);
-template int lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::find_pivot_column_in_row(unsigned int) const;
-template void lp::square_dense_submatrix<lp::mpq, lp::numeric_pair<lp::mpq> >::pivot(unsigned int, lp::lp_settings&);
-#ifdef Z3DEBUG
-template double lp::square_dense_submatrix<double, double>::get_elem(unsigned int, unsigned int) const;
-#endif
-template void lp::square_dense_submatrix<double, double>::apply_from_right(vector<double>&);
-   
-template void  lp::square_dense_submatrix<double, double>::apply_from_left_local<double>(lp::indexed_vector<double>&, lp::lp_settings&);
-template void  lp::square_dense_submatrix<double, double>::apply_from_left_to_vector<double>(vector<double>&);
-template lp::square_dense_submatrix<lp::mpq, lp::mpq>::square_dense_submatrix(lp::square_sparse_matrix<lp::mpq, lp::mpq>*, unsigned int);
-template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::update_parent_matrix(lp::lp_settings&);
-template bool lp::square_dense_submatrix<lp::mpq, lp::mpq>::is_L_matrix() const;
-template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::conjugate_by_permutation(lp::permutation_matrix<lp::mpq, lp::mpq>&);
-template int lp::square_dense_submatrix<lp::mpq, lp::mpq>::find_pivot_column_in_row(unsigned int) const;
-template void lp::square_dense_submatrix<lp::mpq, lp::mpq>::pivot(unsigned int, lp::lp_settings&);
diff --git a/src/math/lp/square_dense_submatrix.h b/src/math/lp/square_dense_submatrix.h
deleted file mode 100644
index 308f9eadc63..00000000000
--- a/src/math/lp/square_dense_submatrix.h
+++ /dev/null
@@ -1,225 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-
-#pragma once
-#include "util/vector.h"
-#include "math/lp/permutation_matrix.h"
-#include <unordered_map>
-#include "math/lp/static_matrix.h"
-#include <set>
-#include <utility>
-#include <string>
-#include <algorithm>
-#include <queue>
-#include "math/lp/indexed_value.h"
-#include "math/lp/indexed_vector.h"
-#include <functional>
-#include "math/lp/lp_settings.h"
-#include "math/lp/eta_matrix.h"
-#include "math/lp/binary_heap_upair_queue.h"
-#include "math/lp/square_sparse_matrix.h"
-namespace lp {
-template <typename T, typename X>
-class square_dense_submatrix : public tail_matrix<T, X> {
-    // the submatrix uses the permutations of the parent matrix to access the elements
-    struct ref {
-        unsigned m_i_offset;
-        square_dense_submatrix & m_s;
-        ref(unsigned i, square_dense_submatrix & s) :
-            m_i_offset((i - s.m_index_start) * s.m_dim), m_s(s){}
-        T & operator[] (unsigned j) {
-            lp_assert(j >= m_s.m_index_start);
-            return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
-        }
-        const T & operator[] (unsigned j) const {
-            lp_assert(j >= m_s.m_index_start);
-            return m_s.m_v[m_i_offset + m_s.adjust_column(j) - m_s.m_index_start];
-        }
-    };
-public:
-    unsigned m_index_start;
-    unsigned m_dim;
-    vector<T> m_v;
-    square_sparse_matrix<T, X> * m_parent;
-    permutation_matrix<T, X>  m_row_permutation;
-    indexed_vector<T> m_work_vector;
-public:
-    permutation_matrix<T, X>  m_column_permutation;
-    bool is_active() const { return m_parent != nullptr; }
-
-    square_dense_submatrix() {}
-
-    square_dense_submatrix (square_sparse_matrix<T, X> *parent_matrix, unsigned index_start);
-
-    void init(square_sparse_matrix<T, X> *parent_matrix, unsigned index_start);
-
-    bool is_dense() const override { return true; }
-    
-    ref operator[] (unsigned i) {
-        lp_assert(i >= m_index_start);
-        lp_assert(i < m_parent->dimension());
-        return ref(i, *this);
-    }
-
-    int find_pivot_column_in_row(unsigned i) const;
-
-    void swap_columns(unsigned i, unsigned j) {
-        if (i != j)
-            m_column_permutation.transpose_from_left(i, j);
-    }
-
-    unsigned adjust_column(unsigned  col)  const{
-        if (col >= m_column_permutation.size())
-            return col;
-        return m_column_permutation.apply_reverse(col);
-    }
-
-    unsigned adjust_column_inverse(unsigned  col)  const{
-        if (col >= m_column_permutation.size())
-            return col;
-        return m_column_permutation[col];
-    }
-    unsigned adjust_row(unsigned row)  const{
-        if (row >= m_row_permutation.size())
-            return row;
-        return m_row_permutation[row];
-    }
-
-    unsigned adjust_row_inverse(unsigned row)  const{
-        if (row >= m_row_permutation.size())
-            return row;
-        return m_row_permutation.apply_reverse(row);
-    }
-
-    void pivot(unsigned i, lp_settings & settings);
-
-    void pivot_row_to_row(unsigned i, unsigned row, lp_settings & settings);;
-
-    void divide_row_by_pivot(unsigned i);
-
-    void update_parent_matrix(lp_settings & settings);
-
-    void update_existing_or_delete_in_parent_matrix_for_row(unsigned i, lp_settings & settings);
-
-    void push_new_elements_to_parent_matrix(lp_settings & settings);
-
-    template <typename L>
-    L row_by_vector_product(unsigned i, const vector<L> & v);
-
-    template <typename L>
-    L column_by_vector_product(unsigned j, const vector<L> & v);
-
-    template <typename L>
-    L row_by_indexed_vector_product(unsigned i, const indexed_vector<L> & v);
-
-    template <typename L>
-    void apply_from_left_local(indexed_vector<L> & w, lp_settings & settings);
-
-    template <typename L>
-    void apply_from_left_to_vector(vector<L> & w);
-
-    bool is_L_matrix() const;
-
-    void apply_from_left_to_T(indexed_vector<T> & w, lp_settings & settings) override {
-        apply_from_left_local(w, settings);
-    }
-
-    
-    
-    void apply_from_right(indexed_vector<T> & w) override {
-#if 1==0
-        indexed_vector<T> wcopy = w;
-        apply_from_right(wcopy.m_data);
-        wcopy.m_index.clear();
-        if (numeric_traits<T>::precise()) {
-            for (unsigned i = 0; i < m_parent->dimension(); i++) {
-                if (!is_zero(wcopy.m_data[i]))
-                    wcopy.m_index.push_back(i);
-            }
-        } else  {
-            for (unsigned i = 0; i < m_parent->dimension(); i++) {
-                T & v = wcopy.m_data[i];
-                if (!lp_settings::is_eps_small_general(v, 1e-14)){
-                    wcopy.m_index.push_back(i);
-                } else {
-                    v = zero_of_type<T>();
-                }
-            }
-        }
-        lp_assert(wcopy.is_OK());
-        apply_from_right(w.m_data);
-        w.m_index.clear();
-        if (numeric_traits<T>::precise()) {
-            for (unsigned i = 0; i < m_parent->dimension(); i++) {
-                if (!is_zero(w.m_data[i]))
-                    w.m_index.push_back(i);
-            }
-        } else  {
-            for (unsigned i = 0; i < m_parent->dimension(); i++) {
-                T & v = w.m_data[i];
-                if (!lp_settings::is_eps_small_general(v, 1e-14)){
-                    w.m_index.push_back(i);
-                } else {
-                    v = zero_of_type<T>();
-                }
-            }
-        }
-#else
-        lp_assert(w.is_OK());
-        lp_assert(m_work_vector.is_OK());
-        m_work_vector.resize(w.data_size());
-        m_work_vector.clear();
-        lp_assert(m_work_vector.is_OK());
-        unsigned end = m_index_start + m_dim;
-        for (unsigned k : w.m_index) {
-            // find j such that k = adjust_row_inverse(j)
-            unsigned j = adjust_row(k);
-            if (j < m_index_start || j >= end) {
-                m_work_vector.set_value(w[k], adjust_column_inverse(j));
-            }  else { // j >= m_index_start and j < end
-                unsigned offset = (j - m_index_start) * m_dim; // this is the row start
-                const T& wv = w[k];
-                for (unsigned col = m_index_start; col < end; col++, offset ++) {
-                    unsigned adj_col =  adjust_column_inverse(col);
-                    m_work_vector.add_value_at_index(adj_col, m_v[offset] * wv);
-                }
-            }
-        }
-        m_work_vector.clean_up();
-        lp_assert(m_work_vector.is_OK());
-        w = m_work_vector;
-#endif
-    }
-    void apply_from_left(vector<X> & w, lp_settings & /*settings*/) override {
-        apply_from_left_to_vector(w);// , settings);
-    }
-
-    void apply_from_right(vector<T> & w) override;
-
-#ifdef Z3DEBUG
-    T get_elem (unsigned i, unsigned j) const override;
-    unsigned row_count() const override { return m_parent->row_count();}
-    unsigned column_count() const override { return row_count();}
-    void set_number_of_rows(unsigned) override {}
-    void set_number_of_columns(unsigned) override {}
-#endif
-    void conjugate_by_permutation(permutation_matrix<T, X> & q);
-};
-}
diff --git a/src/math/lp/square_dense_submatrix_def.h b/src/math/lp/square_dense_submatrix_def.h
deleted file mode 100644
index 3a9006b4d9e..00000000000
--- a/src/math/lp/square_dense_submatrix_def.h
+++ /dev/null
@@ -1,370 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-
-#include "util/vector.h"
-#include "math/lp/square_dense_submatrix.h"
-namespace lp {
-template <typename T, typename X>
-square_dense_submatrix<T, X>::square_dense_submatrix (square_sparse_matrix<T, X> *parent_matrix, unsigned index_start) :
-    m_index_start(index_start),
-    m_dim(parent_matrix->dimension() - index_start),
-    m_v(m_dim * m_dim),
-    m_parent(parent_matrix),
-    m_row_permutation(m_parent->dimension()),
-    m_column_permutation(m_parent->dimension()) {
-    int row_offset = - static_cast<int>(m_index_start);
-    for (unsigned i = index_start; i < parent_matrix->dimension(); i++) {
-        unsigned row = parent_matrix->adjust_row(i);
-        for (auto & iv : parent_matrix->get_row_values(row)) {
-            unsigned j = parent_matrix->adjust_column_inverse(iv.m_index);
-            lp_assert(j>= m_index_start);
-            m_v[row_offset + j] = iv.m_value;
-        }
-        row_offset += m_dim;
-    }
-}
-
-template <typename T, typename X> void square_dense_submatrix<T, X>::init(square_sparse_matrix<T, X> *parent_matrix, unsigned index_start) {
-    m_index_start = index_start;
-    m_dim = parent_matrix->dimension() - index_start;
-    m_v.resize(m_dim * m_dim);
-    m_parent = parent_matrix;
-    m_column_permutation.init(m_parent->dimension());
-    for (unsigned i = index_start; i < parent_matrix->dimension(); i++) {
-        unsigned row = parent_matrix->adjust_row(i);
-        for (auto & iv : parent_matrix->get_row_values(row)) {
-            unsigned j = parent_matrix->adjust_column_inverse(iv.m_index);
-            (*this)[i][j] = iv.m_value;
-        }
-    }
-}
-
-template <typename T, typename X>    int square_dense_submatrix<T, X>::find_pivot_column_in_row(unsigned i) const {
-    int j = -1;
-    T max = zero_of_type<T>();
-    lp_assert(i >= m_index_start);
-    unsigned row_start = (i - m_index_start) * m_dim;
-    for (unsigned k = i; k < m_parent->dimension(); k++) {
-        unsigned col = adjust_column(k); // this is where the column is in the row
-        unsigned offs = row_start + col - m_index_start;
-        T t = abs(m_v[offs]);
-        if (t > max) {
-            j = k;
-            max = t;
-        }
-    }
-    return j;
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::pivot(unsigned i, lp_settings & settings) {
-    divide_row_by_pivot(i);
-    for (unsigned k = i + 1; k < m_parent->dimension(); k++)
-        pivot_row_to_row(i, k, settings);
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::pivot_row_to_row(unsigned i, unsigned row, lp_settings & settings) {
-    lp_assert(i < row);
-    unsigned pj = adjust_column(i); // the pivot column
-    unsigned pjd = pj - m_index_start;
-    unsigned pivot_row_offset = (i-m_index_start)*m_dim;
-    T pivot = m_v[pivot_row_offset + pjd];
-    unsigned row_offset= (row-m_index_start)*m_dim;
-    T m = m_v[row_offset + pjd];
-    lp_assert(!is_zero(pivot));
-    m_v[row_offset + pjd] = -m * pivot; // creating L matrix
-    for (unsigned j = m_index_start; j < m_parent->dimension(); j++) {
-        if (j == pj) {
-            pivot_row_offset++;
-            row_offset++;
-            continue;
-        }
-        auto t = m_v[row_offset] - m_v[pivot_row_offset] * m;
-        if (settings.abs_val_is_smaller_than_drop_tolerance(t)) {
-            m_v[row_offset] = zero_of_type<T>();
-        } else {
-            m_v[row_offset] = t;
-        }
-        row_offset++; pivot_row_offset++;
-        // at the same time we pivot the L too
-    }
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::divide_row_by_pivot(unsigned i) {
-    unsigned pj = adjust_column(i); // the pivot column
-    unsigned irow_offset = (i - m_index_start) * m_dim;
-    T pivot = m_v[irow_offset + pj - m_index_start];
-    lp_assert(!is_zero(pivot));
-    for (unsigned k = m_index_start; k < m_parent->dimension(); k++) {
-        if (k == pj){
-            m_v[irow_offset++] = one_of_type<T>() / pivot; // creating the L matrix diagonal
-            continue;
-        }
-        m_v[irow_offset++] /= pivot;
-    }
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::update_parent_matrix(lp_settings & settings) {
-    for (unsigned i = m_index_start; i < m_parent->dimension(); i++)
-        update_existing_or_delete_in_parent_matrix_for_row(i, settings);
-    push_new_elements_to_parent_matrix(settings);
-    for (unsigned i = m_index_start; i < m_parent->dimension(); i++)
-        m_parent->set_max_in_row(m_parent->adjust_row(i));
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::update_existing_or_delete_in_parent_matrix_for_row(unsigned i, lp_settings & settings) {
-    bool diag_updated = false;
-    unsigned ai = m_parent->adjust_row(i);
-    auto & row_vals = m_parent->get_row_values(ai);
-    for (unsigned k = 0; k < row_vals.size(); k++) {
-        auto & iv = row_vals[k];
-        unsigned j = m_parent->adjust_column_inverse(iv.m_index);
-        if (j < i) {
-            m_parent->remove_element(row_vals, iv);
-            k--;
-        } else if (i == j) {
-            m_parent->m_columns[iv.m_index].m_values[iv.m_other].set_value(iv.m_value = one_of_type<T>());
-            diag_updated = true;
-        } else { // j > i
-            T & v = (*this)[i][j];
-            if (settings.abs_val_is_smaller_than_drop_tolerance(v)) {
-                m_parent->remove_element(row_vals, iv);
-                k--;
-            } else {
-                m_parent->m_columns[iv.m_index].m_values[iv.m_other].set_value(iv.m_value = v);
-                v = zero_of_type<T>(); // only new elements are left above the diagonal
-            }
-        }
-    }
-    if (!diag_updated) {
-        unsigned aj = m_parent->adjust_column(i);
-        m_parent->add_new_element(ai, aj, one_of_type<T>());
-    }
-}
-
-template <typename T, typename X>    void square_dense_submatrix<T, X>::push_new_elements_to_parent_matrix(lp_settings & settings) {
-    for (unsigned i = m_index_start; i < m_parent->dimension() - 1; i++) {
-        unsigned ai = m_parent->adjust_row(i);
-        for (unsigned j = i + 1; j < m_parent->dimension(); j++) {
-            T & v = (*this)[i][j];
-            if (!settings.abs_val_is_smaller_than_drop_tolerance(v)) {
-                unsigned aj = m_parent->adjust_column(j);
-                m_parent->add_new_element(ai, aj, v);
-            }
-            v = zero_of_type<T>(); // leave only L elements now
-        }
-    }
-}
-template <typename T, typename X>
-template <typename L>
-L square_dense_submatrix<T, X>::row_by_vector_product(unsigned i, const vector<L> & v) {
-    lp_assert(i >= m_index_start);
-
-    unsigned row_in_subm = i - m_index_start;
-    unsigned row_offset = row_in_subm * m_dim;
-    L r = zero_of_type<L>();
-    for (unsigned j = 0; j < m_dim; j++)
-        r += m_v[row_offset + j] * v[adjust_column_inverse(m_index_start + j)];
-    return r;
-}
-
-template <typename T, typename X>
-template <typename L>
-L square_dense_submatrix<T, X>::column_by_vector_product(unsigned j, const vector<L> & v) {
-    lp_assert(j >= m_index_start);
-
-    unsigned offset = j - m_index_start;
-    L r = zero_of_type<L>();
-    for (unsigned i = 0; i < m_dim; i++, offset += m_dim)
-        r += m_v[offset] * v[adjust_row_inverse(m_index_start + i)];
-    return r;
-}
-template <typename T, typename X>
-template <typename L>
-L square_dense_submatrix<T, X>::row_by_indexed_vector_product(unsigned i, const indexed_vector<L> & v) {
-    lp_assert(i >= m_index_start);
-
-    unsigned row_in_subm = i - m_index_start;
-    unsigned row_offset = row_in_subm * m_dim;
-    L r = zero_of_type<L>();
-    for (unsigned j = 0; j < m_dim; j++)
-        r += m_v[row_offset + j] * v[adjust_column_inverse(m_index_start + j)];
-    return r;
-}
-template <typename T, typename X>
-template <typename L>
-void square_dense_submatrix<T, X>::apply_from_left_local(indexed_vector<L> & w, lp_settings & settings) {
-#ifdef Z3DEBUG
-    // dense_matrix<T, X> deb(*this);
-    // vector<L>  deb_w(w.m_data.size());
-    // for (unsigned i = 0; i < w.m_data.size(); i++)
-    //     deb_w[i] = w[i];
-
-    // deb.apply_from_left(deb_w);
-#endif // use indexed vector here
-
-#ifndef DO_NOT_USE_INDEX
-    vector<L> t(m_parent->dimension(), zero_of_type<L>());
-    for (auto k : w.m_index) {
-        unsigned j = adjust_column(k); // k-th element will contribute only to column j
-        if (j < m_index_start || j >= this->m_index_start +  this->m_dim) { // it is a unit matrix outside 
-            t[adjust_row_inverse(j)] = w[k];
-        } else {
-            const L & v = w[k];
-            for (unsigned i = 0; i < m_dim; i++) {
-                unsigned row = adjust_row_inverse(m_index_start + i);
-                unsigned offs = i * m_dim + j - m_index_start;
-                t[row] += m_v[offs] * v;
-            }
-        }
-    }
-    w.m_index.clear();
-    for (unsigned i = 0; i < m_parent->dimension(); i++) {
-        const L & v = t[i];
-        if (!settings.abs_val_is_smaller_than_drop_tolerance(v)){
-            w.m_index.push_back(i);
-            w.m_data[i] = v;
-        } else {
-            w.m_data[i] = zero_of_type<L>();
-        }
-    }
-#else
-    vector<L> t(m_parent->dimension());
-    for (unsigned i = 0; i < m_index_start; i++) {
-        t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
-    }
-    for (unsigned i = m_index_start; i < m_parent->dimension(); i++){
-        t[adjust_row_inverse(i)] = row_by_indexed_vector_product(i, w);
-    }
-    for (unsigned i = 0; i < m_parent->dimension(); i++) {
-        w.set_value(t[i], i);
-    }
-    for (unsigned i = 0; i < m_parent->dimension(); i++) {
-        const L & v = t[i];
-        if (!is_zero(v))
-            w.m_index.push_back(i);
-        w.m_data[i] = v;
-    }
-#endif
-#ifdef Z3DEBUG
-    // cout << "w final" << endl;
-    // print_vector(w.m_data);
-    //        lp_assert(vectors_are_equal<T>(deb_w, w.m_data));
-    // lp_assert(w.is_OK());
-#endif
-}
-
-template <typename T, typename X>
-template <typename L>
-void square_dense_submatrix<T, X>::apply_from_left_to_vector(vector<L> & w) {
-    // lp_settings & settings) {
-    // dense_matrix<T, L> deb(*this);
-    // vector<L>  deb_w(w);
-    // deb.apply_from_left_to_X(deb_w, settings);
-    // // cout << "deb" << endl;
-    // // print_matrix(deb);
-    // // cout << "w" << endl;
-    // // print_vector(w.m_data);
-    // // cout << "deb_w" << endl;
-    // // print_vector(deb_w);
-    vector<L> t(m_parent->dimension());
-    for (unsigned i = 0; i < m_index_start; i++) {
-        t[adjust_row_inverse(i)] = w[adjust_column_inverse(i)];
-    }
-    for (unsigned i = m_index_start; i < m_parent->dimension(); i++){
-        t[adjust_row_inverse(i)] = row_by_vector_product(i, w);
-    }
-    for (unsigned i = 0; i < m_parent->dimension(); i++) {
-        w[i] = t[i];
-    }
-#ifdef Z3DEBUG
-    // cout << "w final" << endl;
-    // print_vector(w.m_data);
-    //  lp_assert(vectors_are_equal<L>(deb_w, w));
-#endif
-}
-
-template <typename T, typename X>    bool square_dense_submatrix<T, X>::is_L_matrix() const {
-#ifdef Z3DEBUG
-    lp_assert(m_row_permutation.is_identity());
-    for (unsigned i = 0; i < m_parent->dimension(); i++) {
-        if (i < m_index_start) {
-            lp_assert(m_column_permutation[i] == i);
-            continue;
-        }
-        unsigned row_offs = (i-m_index_start)*m_dim;
-        for (unsigned k = 0; k < m_dim; k++) {
-            unsigned j = m_index_start + k;
-            unsigned jex = adjust_column_inverse(j);
-            if (jex > i) {
-                lp_assert(is_zero(m_v[row_offs + k]));
-            } else if (jex == i) {
-                lp_assert(!is_zero(m_v[row_offs + k]));
-            }
-        }
-    }
-#endif
-    return true;
-}
-
-template <typename T, typename X> void square_dense_submatrix<T, X>::apply_from_right(vector<T> & w) {
-#ifdef Z3DEBUG
-    // dense_matrix<T, X> deb(*this);
-    // vector<T>  deb_w(w);
-    // deb.apply_from_right(deb_w);
-#endif
-    vector<T> t(w.size());
-
-    for (unsigned j = 0; j < m_index_start; j++) {
-        t[adjust_column_inverse(j)] = w[adjust_row_inverse(j)];
-    }
-    unsigned end = m_index_start + m_dim;
-    for (unsigned j = end; j < m_parent->dimension(); j++) {
-        t[adjust_column_inverse(j)] = w[adjust_row_inverse(j)];
-    }
-    for (unsigned j = m_index_start; j < end; j++) {
-        t[adjust_column_inverse(j)] = column_by_vector_product(j, w);
-    }
-    w = t;
-#ifdef Z3DEBUG
-    //  lp_assert(vector_are_equal<T>(deb_w, w));
-#endif
-}
-
-
-
-
-#ifdef Z3DEBUG
-
-template <typename T, typename X>    T square_dense_submatrix<T, X>::get_elem (unsigned i, unsigned j) const {
-    i = adjust_row(i);
-    j = adjust_column(j);
-    if (i < m_index_start || j < m_index_start)
-        return i == j? one_of_type<T>() : zero_of_type<T>();
-    unsigned offs = (i - m_index_start)* m_dim + j - m_index_start;
-    return  m_v[offs];
-}
-
-#endif
-template <typename T, typename X> void square_dense_submatrix<T, X>::conjugate_by_permutation(permutation_matrix<T, X> & q) {
-    m_row_permutation.multiply_by_permutation_from_left(q);
-    m_column_permutation.multiply_by_reverse_from_right(q);
-}
-}
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index c386cf002d9..2e8faeaedc5 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -48,7 +48,6 @@
 #include "test/lp/gomory_test.h"
 #include "math/lp/matrix.h"
 #include "math/lp/hnf.h"
-#include "math/lp/square_sparse_matrix_def.h"
 #include "math/lp/general_matrix.h"
 #include "math/lp/lp_bound_propagator.h"
 #include "math/lp/nla_solver.h"
@@ -56,6 +55,8 @@
 #include "math/lp/cross_nested.h"
 #include "math/lp/int_cube.h"
 #include "math/lp/emonics.h"
+#include "math/lp/static_matrix.h"
+#include "math/lp/dense_matrix.h"
 
 bool my_white_space(const char & a) {
     return a == ' ' || a == '\t';
@@ -388,122 +389,6 @@ struct simple_column_namer:public column_namer
 };
 
 
-template <typename T, typename X>
-void test_matrix(square_sparse_matrix<T, X> & a) {
-    auto m = a.dimension();
-
-    // copy a to b in the reversed order
-    square_sparse_matrix<T, X> b(m, m);
-    std::cout << "copy b to a"<< std::endl;
-    for (int row = m - 1; row >= 0; row--)
-        for (int col = m - 1; col >= 0; col --) {
-            b(row, col) = (T const&) a(row, col);
-        }
-
-
-    std::cout << "zeroing b in the reverse order"<< std::endl;
-    for (int row = m - 1; row >= 0; row--)
-        for (int col = m - 1; col >= 0; col --)
-            b.set(row, col, T(0));
-
-
-
-    for (unsigned row = 0; row < m; row ++)
-        for (unsigned col = 0; col < m; col ++)
-            a.set(row, col, T(0));
-
-
-    unsigned i = my_random() % m;
-    unsigned j = my_random() % m;
-
-    auto t = T(1);
-
-    a.set(i, j, t);
-
-    lp_assert(a.get(i, j) == t);
-
-    unsigned j1;
-    if (j < m - 1) {
-        j1 = m - 1;
-        a.set(i, j1, T(2));
-    }
-}
-
-void tst1() {
-    std::cout << "testing the minimal matrix with 1 row and 1 column" << std::endl;
-    square_sparse_matrix<double, double> m0(1, 1);
-    m0.set(0, 0, 1);
-    // print_matrix(m0);
-    m0.set(0, 0, 0);
-    // print_matrix(m0);
-    test_matrix(m0);
-
-    unsigned rows = 2;
-    square_sparse_matrix<double, double> m(rows, rows);
-    std::cout << "setting m(0,1)=" << std::endl;
-
-    m.set(0, 1,  11);
-    m.set(0, 0,  12);
-
-    // print_matrix<double>(m);
-
-    test_matrix(m);
-
-    square_sparse_matrix<double, double> m1(2, 2);
-    m1.set(0, 0, 2);
-    m1.set(1, 0, 3);
-    // print_matrix(m1);
-    std::cout << " zeroing matrix 2 by 2" << std::endl;
-    m1.set(0, 0, 0);
-    m1.set(1, 0, 0);
-    // print_matrix(m1);
-
-    test_matrix(m1);
-
-
-    std::cout << "printing zero matrix 3 by 1" << std::endl;
-    square_sparse_matrix<double, double> m2(3, 3);
-    // print_matrix(m2);
-
-    m2.set(0, 0, 1);
-    m2.set(2, 0, 2);
-    std::cout << "printing  matrix 3 by 1 with a gap" << std::endl;
-    // print_matrix(m2);
-
-    test_matrix(m2);
-
-    square_sparse_matrix<double, double> m10by9(10, 10);
-    m10by9.set(0, 1, 1);
-
-    m10by9(0, 1) = 4;
-
-    double test = m10by9(0, 1);
-
-    std::cout << "got " << test << std::endl;
-
-
-    m10by9.set(0, 8, 8);
-    m10by9.set(3, 4, 7);
-    m10by9.set(3, 2, 5);
-    m10by9.set(3, 8, 99);
-    m10by9.set(3, 2, 6);
-    m10by9.set(1, 8, 9);
-    m10by9.set(4, 0, 40);
-    m10by9.set(0, 0, 10);
-
-    std::cout << "printing  matrix 10 by 9" << std::endl;
-    // print_matrix(m10by9);
-
-
-    test_matrix(m10by9);
-    std::cout <<"zeroing m10by9\n";
-#ifdef Z3DEBUG
-    for (unsigned int i = 0; i < m10by9.dimension(); i++)
-        for (unsigned int j = 0; j < m10by9.column_count(); j++)
-            m10by9.set(i, j, 0);
-#endif
-    // print_matrix(m10by9);
-}
 
 vector<int> allocate_basis_heading(unsigned count) { // the rest of initialization will be handled by lu_QR
     vector<int> basis_heading(count, -1);
@@ -555,12 +440,6 @@ void test_small_lu(lp_settings & settings) {
 
 #endif
 
-void fill_long_row(square_sparse_matrix<double, double> &m, int i) {
-    int n = m.dimension();
-    for (int j = 0; j < n; j ++) {
-        m (i, (j + i) % n) = j * j;
-    }
-}
 
 void fill_long_row(static_matrix<double, double> &m, int i) {
     int n = m.column_count();
@@ -570,13 +449,6 @@ void fill_long_row(static_matrix<double, double> &m, int i) {
 }
 
 
-void fill_long_row_exp(square_sparse_matrix<double, double> &m, int i) {
-    int n = m.dimension();
-
-    for (int j = 0; j < n; j ++) {
-        m(i, j) = my_random() % 20;
-    }
-}
 
 void fill_long_row_exp(static_matrix<double, double> &m, int i) {
     int n = m.column_count();
@@ -586,26 +458,9 @@ void fill_long_row_exp(static_matrix<double, double> &m, int i) {
     }
 }
 
-void fill_larger_square_sparse_matrix_exp(square_sparse_matrix<double, double> & m){
-    for ( unsigned i = 0; i < m.dimension(); i++ )
-        fill_long_row_exp(m, i);
-}
-
-void fill_larger_square_sparse_matrix_exp(static_matrix<double, double> & m){
-    for ( unsigned i = 0; i < m.row_count(); i++ )
-        fill_long_row_exp(m, i);
-}
 
 
-void fill_larger_square_sparse_matrix(square_sparse_matrix<double, double> & m){
-    for ( unsigned i = 0; i < m.dimension(); i++ )
-        fill_long_row(m, i);
-}
 
-void fill_larger_square_sparse_matrix(static_matrix<double, double> & m){
-    for ( unsigned i = 0; i < m.row_count(); i++ )
-        fill_long_row(m, i);
-}
 
 
 int perm_id = 0;
@@ -616,11 +471,6 @@ int perm_id = 0;
 
 
 
-void init_b(vector<double> & b, square_sparse_matrix<double, double> & m, vector<double>& x) {
-    for (unsigned i = 0; i < m.dimension(); i++) {
-        b.push_back(m.dot_product_with_row(i, x));
-    }
-}
 
 void init_b(vector<double> & b, static_matrix<double, double> & m, vector<double> & x) {
     for (unsigned i = 0; i < m.row_count(); i++) {
@@ -671,7 +521,7 @@ void test_lp_1() {
     m(2, 0) =  2;  m(2, 1) = -1; m(2, 2) = 2;  m(2, 5) = 1;
     m(3, 0) =  2;  m(3, 1) =  3; m(3, 2) = -1; m(3, 6) = 1;
 #ifdef Z3DEBUG
-    print_matrix(m, std::cout);
+    //print_matrix(m, std::cout);
 #endif
     vector<double> x_star(7);
     x_star[0] = 0; x_star[1] = 0; x_star[2] = 0;
@@ -724,213 +574,9 @@ void test_lp_primal_core_solver() {
 }
 
 
-#ifdef Z3DEBUG
-template <typename T, typename X>
-void test_swap_rows_with_permutation(square_sparse_matrix<T, X>& m){
-    std::cout << "testing swaps" << std::endl;
-    unsigned dim = m.row_count();
-    dense_matrix<double, double> original(&m);
-    permutation_matrix<double, double> q(dim);
-    print_matrix(m, std::cout);
-    lp_assert(original == q * m);
-    for (int i = 0; i < 100; i++) {
-        unsigned row1 = my_random() % dim;
-        unsigned row2 = my_random() % dim;
-        if (row1 == row2) continue;
-        std::cout << "swap " << row1 << " " << row2 << std::endl;
-        m.swap_rows(row1, row2);
-        q.transpose_from_left(row1, row2);
-        lp_assert(original == q * m);
-        print_matrix(m, std::cout);
-        std::cout << std::endl;
-    }
-}
-#endif
-template <typename T, typename X>
-void fill_matrix(square_sparse_matrix<T, X>& m); // forward definition
-#ifdef Z3DEBUG
-template <typename T, typename X>
-void test_swap_cols_with_permutation(square_sparse_matrix<T, X>& m){
-    std::cout << "testing swaps" << std::endl;
-    unsigned dim = m.row_count();
-    dense_matrix<double, double> original(&m);
-    permutation_matrix<double, double> q(dim);
-    print_matrix(m, std::cout);
-    lp_assert(original == q * m);
-    for (int i = 0; i < 100; i++) {
-        unsigned row1 = my_random() % dim;
-        unsigned row2 = my_random() % dim;
-        if (row1 == row2) continue;
-        std::cout << "swap " << row1 << " " << row2 << std::endl;
-        m.swap_rows(row1, row2);
-        q.transpose_from_right(row1, row2);
-        lp_assert(original == q * m);
-        print_matrix(m, std::cout);
-        std::cout << std::endl;
-    }
-}
-
-
-template <typename T, typename X>
-void test_swap_rows(square_sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
-    std::cout << "test_swap_rows(" << i0 << "," << i1 << ")" << std::endl;
-    square_sparse_matrix<T, X> mcopy(m.dimension(), 0);
-    for (unsigned i = 0; i  < m.dimension(); i++)
-        for (unsigned j = 0; j < m.dimension(); j++) {
-            mcopy(i, j)= m(i, j);
-        }
-    std::cout << "swapping rows "<< i0 << "," << i1 << std::endl;
-    m.swap_rows(i0, i1);
-
-    for (unsigned j = 0; j < m.dimension(); j++) {
-        lp_assert(mcopy(i0, j) == m(i1, j));
-        lp_assert(mcopy(i1, j) == m(i0, j));
-    }
-}
-template <typename T, typename X>
-void test_swap_columns(square_sparse_matrix<T, X>& m, unsigned i0, unsigned i1){
-    std::cout << "test_swap_columns(" << i0 << "," << i1 << ")" << std::endl;
-    square_sparse_matrix<T, X> mcopy(m.dimension(), 0); // the second argument does not matter
-    for (unsigned i = 0; i  < m.dimension(); i++)
-        for (unsigned j = 0; j < m.dimension(); j++) {
-            mcopy(i, j)= m(i, j);
-        }
-    m.swap_columns(i0, i1);
-
-    for (unsigned j = 0; j < m.dimension(); j++) {
-        lp_assert(mcopy(j, i0) == m(j, i1));
-        lp_assert(mcopy(j, i1) == m(j, i0));
-    }
-
-    for (unsigned i = 0; i  < m.dimension(); i++) {
-        if (i == i0 || i == i1)
-            continue;
-        for (unsigned j = 0; j < m.dimension(); j++) {
-            lp_assert(mcopy(j, i)== m(j, i));
-        }
-    }
-}
-#endif
-
-template <typename T, typename X>
-void fill_matrix(square_sparse_matrix<T, X>& m){
-    int v = 0;
-    for (int i = m.dimension() - 1; i >= 0; i--) {
-        for (int j = m.dimension() - 1; j >=0; j--){
-            m(i, j) = v++;
-        }
-    }
-}
-
-void test_pivot_like_swaps_and_pivot(){
-    square_sparse_matrix<double, double> m(10, 10);
-    fill_matrix(m);
-    // print_matrix(m);
-    // pivot at 2,7
-    m.swap_columns(0, 7);
-    // print_matrix(m);
-    m.swap_rows(2, 0);
-    // print_matrix(m);
-    for (unsigned i = 1; i < m.dimension(); i++) {
-        m(i, 0) = 0;
-    }
-    // print_matrix(m);
-
-    // say pivot at 3,4
-    m.swap_columns(1, 4);
-    // print_matrix(m);
-    m.swap_rows(1, 3);
-    // print_matrix(m);
-
-    vector<double> row;
-    double alpha = 2.33;
-    unsigned pivot_row = 1;
-    unsigned target_row = 2;
-    unsigned pivot_row_0 = 3;
-    double beta = 3.1;
-    m(target_row, 3) = 0;
-    m(target_row, 5) = 0;
-    m(pivot_row, 6) = 0;
-#ifdef Z3DEBUG
-    print_matrix(m, std::cout);
-#endif
-
-    for (unsigned j = 0; j < m.dimension(); j++) {
-        row.push_back(m(target_row, j) + alpha * m(pivot_row, j) + beta * m(pivot_row_0, j));
-    }
-
-    for (auto & t : row) {
-        std::cout << t << ",";
-    }
-    std::cout << std::endl;
-    lp_settings settings;
-    m.pivot_row_to_row(pivot_row, alpha, target_row, settings);
-    m.pivot_row_to_row(pivot_row_0, beta, target_row, settings);
-    //  print_matrix(m);
-    for (unsigned j = 0; j < m.dimension(); j++) {
-        lp_assert(abs(row[j] - m(target_row, j)) < 0.00000001);
-    }
-}
-
-#ifdef Z3DEBUG
-void test_swap_rows() {
-    square_sparse_matrix<double, double> m(10, 10);
-    fill_matrix(m);
-    // print_matrix(m);
-    test_swap_rows(m, 3, 5);
-
-    test_swap_rows(m, 1, 3);
-
-
-    test_swap_rows(m, 1, 3);
-
-    test_swap_rows(m, 1, 7);
-
-    test_swap_rows(m, 3, 7);
 
-    test_swap_rows(m, 0, 7);
 
-    m(0, 4) = 1;
-    // print_matrix(m);
-    test_swap_rows(m, 0, 7);
 
-    // go over some corner cases
-    square_sparse_matrix<double, double> m0(2, 2);
-    test_swap_rows(m0, 0, 1);
-    m0(0, 0) = 3;
-    test_swap_rows(m0, 0, 1);
-    m0(1, 0) = 3;
-    test_swap_rows(m0, 0, 1);
-
-
-    square_sparse_matrix<double, double> m1(10, 10);
-    test_swap_rows(m1, 0, 1);
-    m1(0, 0) = 3;
-    test_swap_rows(m1, 0, 1);
-    m1(1, 0) = 3;
-    m1(0, 3) = 5;
-    m1(1, 3) = 4;
-    m1(1, 8) = 8;
-    m1(1, 9) = 8;
-
-    test_swap_rows(m1, 0, 1);
-
-    square_sparse_matrix<double, double> m2(3, 3);
-    test_swap_rows(m2, 0, 1);
-    m2(0, 0) = 3;
-    test_swap_rows(m2, 0, 1);
-    m2(2, 0) = 3;
-    test_swap_rows(m2, 0, 2);
-}
-
-void fill_uniformly(square_sparse_matrix<double, double> & m, unsigned dim) {
-    int v = 0;
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            m(i, j) = v++;
-        }
-    }
-}
 
 void fill_uniformly(dense_matrix<double, double> & m, unsigned dim) {
     int v = 0;
@@ -941,146 +587,8 @@ void fill_uniformly(dense_matrix<double, double> & m, unsigned dim) {
     }
 }
 
-void square_sparse_matrix_with_permutations_test() {
-    unsigned dim = 4;
-    square_sparse_matrix<double, double> m(dim, dim);
-    fill_uniformly(m, dim);
-    dense_matrix<double, double> dm(dim, dim);
-    fill_uniformly(dm, dim);
-    dense_matrix<double, double> dm0(dim, dim);
-    fill_uniformly(dm0, dim);
-    permutation_matrix<double, double> q0(dim);
-    q0[0] = 1;
-    q0[1] = 0;
-    q0[2] = 3;
-    q0[3] = 2;
-    permutation_matrix<double, double> q1(dim);
-    q1[0] = 1;
-    q1[1] = 2;
-    q1[2] = 3;
-    q1[3] = 0;
-    permutation_matrix<double, double> p0(dim);
-    p0[0] = 1;
-    p0[1] = 0;
-    p0[2] = 3;
-    p0[3] = 2;
-    permutation_matrix<double, double> p1(dim);
-    p1[0] = 1;
-    p1[1] = 2;
-    p1[2] = 3;
-    p1[3] = 0;
-
-    m.multiply_from_left(q0);
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            lp_assert(m(i, j) == dm0.get_elem(q0[i], j));
-        }
-    }
-
-    auto q0_dm = q0 * dm;
-    lp_assert(m == q0_dm);
-
-    m.multiply_from_left(q1);
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], j));
-        }
-    }
-
-
-    auto q1_q0_dm = q1 * q0_dm;
-
-    lp_assert(m == q1_q0_dm);
-
-    m.multiply_from_right(p0);
-
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p0[j]));
-        }
-    }
-
-    auto q1_q0_dm_p0 = q1_q0_dm * p0;
-
-    lp_assert(m == q1_q0_dm_p0);
-
-    m.multiply_from_right(p1);
-
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p0[j]]));
-        }
-    }
-
-    auto q1_q0_dm_p0_p1 = q1_q0_dm_p0 * p1;
-    lp_assert(m == q1_q0_dm_p0_p1);
-
-    m.multiply_from_right(p1);
-    for (unsigned i = 0; i < dim; i++) {
-        for (unsigned j = 0; j < dim; j++) {
-            lp_assert(m(i, j) == dm0.get_elem(q0[q1[i]], p1[p1[p0[j]]]));
-        }
-    }
-    auto q1_q0_dm_p0_p1_p1 = q1_q0_dm_p0_p1 * p1;
-
-    lp_assert(m == q1_q0_dm_p0_p1_p1);
-}
-
-void test_swap_columns() {
-    square_sparse_matrix<double, double> m(10, 10);
-    fill_matrix(m);
-    // print_matrix(m);
-
-    test_swap_columns(m, 3, 5);
-
-    test_swap_columns(m, 1, 3);
-
-    test_swap_columns(m, 1, 3);
-
-    // print_matrix(m);
-    test_swap_columns(m, 1, 7);
-
-    test_swap_columns(m, 3, 7);
-
-    test_swap_columns(m, 0, 7);
-
-    test_swap_columns(m, 0, 7);
-
-    // go over some corner cases
-    square_sparse_matrix<double, double> m0(2, 2);
-    test_swap_columns(m0, 0, 1);
-    m0(0, 0) = 3;
-    test_swap_columns(m0, 0, 1);
-    m0(0, 1) = 3;
-    test_swap_columns(m0, 0, 1);
-
-
-    square_sparse_matrix<double, double> m1(10, 10);
-    test_swap_columns(m1, 0, 1);
-    m1(0, 0) = 3;
-    test_swap_columns(m1, 0, 1);
-    m1(0, 1) = 3;
-    m1(3, 0) = 5;
-    m1(3, 1) = 4;
-    m1(8, 1) = 8;
-    m1(9, 1) = 8;
-
-    test_swap_columns(m1, 0, 1);
-
-    square_sparse_matrix<double, double> m2(3, 3);
-    test_swap_columns(m2, 0, 1);
-    m2(0, 0) = 3;
-    test_swap_columns(m2, 0, 1);
-    m2(0, 2) = 3;
-    test_swap_columns(m2, 0, 2);
-}
-
 
 
-void test_swap_operations() {
-    test_swap_rows();
-    test_swap_columns();
-}
 
 void test_dense_matrix() {
     dense_matrix<double, double> d(3, 2);
@@ -1110,7 +618,7 @@ void test_dense_matrix() {
     p2.set_elem(1, 0, 1);
     auto c2 = d * p2;
 }
-#endif
+
 
 
 
@@ -1471,59 +979,6 @@ bool contains(std::string const & s, char const * pattern) {
 
 
 
-
-void test_init_U() {
-    static_matrix<double, double> m(3, 7);
-    m(0, 0) = 10; m(0, 1) = 11; m(0, 2) = 12; m(0, 3) = 13; m(0, 4) = 14;
-    m(1, 0) = 20; m(1, 1) = 21; m(1, 2) = 22; m(1, 3) = 23; m(1, 5) = 24;
-    m(2, 0) = 30; m(2, 1) = 31; m(2, 2) = 32; m(2, 3) = 33; m(2, 6) = 34;
-#ifdef Z3DEBUG
-    print_matrix(m, std::cout);
-#endif
-    vector<unsigned> basis(3);
-    basis[0] = 1;
-    basis[1] = 2;
-    basis[2] = 4;
-
-    square_sparse_matrix<double, double> u(m, basis);
-
-    for (unsigned i = 0; i < 3; i++) {
-        for (unsigned j = 0; j < 3; j ++) {
-            lp_assert(m(i, basis[j]) == u(i, j));
-        }
-    }
-
-    // print_matrix(m);
-    // print_matrix(u);
-}
-
-void test_replace_column() {
-    square_sparse_matrix<double, double> m(10, 10);
-    fill_matrix(m);
-    m.swap_columns(0, 7);
-    m.swap_columns(6, 3);
-    m.swap_rows(2, 0);
-    for (unsigned i = 1; i < m.dimension(); i++) {
-        m(i, 0) = 0;
-    }
-
-    indexed_vector<double> w(m.dimension());
-    for (unsigned i = 0; i < m.dimension(); i++) {
-        w.set_value(i % 3, i);
-    }
-
-    lp_settings settings;
-
-    for (unsigned column_to_replace = 0;  column_to_replace < m.dimension(); column_to_replace ++) {
-        m.replace_column(column_to_replace, w, settings);
-        for (unsigned i = 0; i < m.dimension(); i++) {
-            lp_assert(abs(w[i] - m(i, column_to_replace)) < 0.00000001);
-        }
-    }
-}
-
-
-
 void setup_args_parser(argument_parser & parser) {
     parser.add_option_with_help_string("-monics", "test emonics");
     parser.add_option_with_help_string("-nex_order", "test nex order");
@@ -1544,8 +999,6 @@ void setup_args_parser(argument_parser & parser) {
     parser.add_option_with_help_string("-xyz_sample", "run a small interactive scenario");
     parser.add_option_with_after_string_with_help("--density", "the percentage of non-zeroes in the matrix below which it is not dense");
     parser.add_option_with_after_string_with_help("--harris_toler", "harris tolerance");
-    parser.add_option_with_help_string("--test_swaps", "test row swaps with a permutation");
-    parser.add_option_with_help_string("--test_perm", "test permutations");
     parser.add_option_with_after_string_with_help("--checklu", "the file name for lu checking");
     parser.add_option_with_after_string_with_help("--partial_pivot", "the partial pivot constant, a number somewhere between 10 and 100");
     parser.add_option_with_after_string_with_help("--percent_for_enter", "which percent of columns check for entering column");
@@ -1871,38 +1324,6 @@ void check_lu_from_file(std::string lufile_name) {
     lp_assert(false);
 }
 
-void test_square_dense_submatrix() {
-    std::cout << "testing square_dense_submatrix" << std::endl;
-    unsigned parent_dim = 7;
-    square_sparse_matrix<double, double> parent(parent_dim, 0);
-    fill_matrix(parent);
-    unsigned index_start = 3;
-    square_dense_submatrix<double, double> d;
-    d.init(&parent, index_start);
-    for (unsigned i = index_start; i < parent_dim; i++)
-        for (unsigned j = index_start; j < parent_dim; j++)
-            d[i][j] = i*3+j*2;
-#ifdef Z3DEBUG
-    unsigned dim = parent_dim - index_start;
-    dense_matrix<double, double> m(dim, dim);
-    for (unsigned i = index_start; i < parent_dim; i++)
-        for (unsigned j = index_start; j < parent_dim; j++)
-            m[i-index_start][j-index_start] = d[i][j];
-    print_matrix(&m, std::cout);
-#endif
-    for (unsigned i = index_start; i < parent_dim; i++)
-        for (unsigned j = index_start; j < parent_dim; j++)
-            d[i][j] = d[j][i];
-#ifdef Z3DEBUG
-    for (unsigned i = index_start; i < parent_dim; i++)
-        for (unsigned j = index_start; j < parent_dim; j++)
-            m[i-index_start][j-index_start] = d[i][j];
-
-    print_matrix(&m, std::cout);
-    std::cout << std::endl;
-#endif
-}
-
 
 
 void print_st(lp_status status) {
@@ -2958,8 +2379,3 @@ void test_lp_local(int argn, char**argv) {
 void tst_lp(char ** argv, int argc, int& i) {
     lp::test_lp_local(argc - 2, argv + 2);
 }
-#ifdef Z3DEBUG
-namespace lp {
-template void print_matrix<general_matrix>(general_matrix&, std::ostream&);
-}
-#endif