diff --git a/src/math/lp/lp_primal_simplex.cpp b/src/math/lp/lp_primal_simplex.cpp
deleted file mode 100644
index 634f52900..000000000
--- a/src/math/lp/lp_primal_simplex.cpp
+++ /dev/null
@@ -1,35 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#include <utility>
-#include <memory>
-#include <string>
-#include "util/vector.h"
-#include <functional>
-#include "math/lp/lp_primal_simplex_def.h"
-template bool lp::lp_primal_simplex<double, double>::bounds_hold(std::unordered_map<std::string, double, std::hash<std::string>, std::equal_to<std::string>, std::allocator<std::pair<std::string const, double> > > const&);
-template bool lp::lp_primal_simplex<double, double>::row_constraints_hold(std::unordered_map<std::string, double, std::hash<std::string>, std::equal_to<std::string>, std::allocator<std::pair<std::string const, double> > > const&);
-template double lp::lp_primal_simplex<double, double>::get_current_cost() const;
-template double lp::lp_primal_simplex<double, double>::get_column_value(unsigned int) const;
-template lp::lp_primal_simplex<double, double>::~lp_primal_simplex();
-template lp::lp_primal_simplex<lp::mpq, lp::mpq>::~lp_primal_simplex();
-template lp::mpq lp::lp_primal_simplex<lp::mpq, lp::mpq>::get_current_cost() const;
-template lp::mpq lp::lp_primal_simplex<lp::mpq, lp::mpq>::get_column_value(unsigned int) const;
-template void lp::lp_primal_simplex<double, double>::find_maximal_solution();
-template void lp::lp_primal_simplex<lp::mpq, lp::mpq>::find_maximal_solution();
diff --git a/src/math/lp/lp_primal_simplex.h b/src/math/lp/lp_primal_simplex.h
deleted file mode 100644
index 77e12d088..000000000
--- a/src/math/lp/lp_primal_simplex.h
+++ /dev/null
@@ -1,106 +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 <unordered_map>
-#include <string>
-#include <algorithm>
-#include "math/lp/lp_utils.h"
-#include "math/lp/column_info.h"
-#include "math/lp/lp_primal_core_solver.h"
-#include "math/lp/lp_solver.h"
-namespace lp {
-template <typename T, typename X>
-class lp_primal_simplex: public lp_solver<T, X> {
-    lp_primal_core_solver<T, X> * m_core_solver;
-    vector<X> m_lower_bounds;
-private:
-    unsigned original_rows() { return this->m_external_rows_to_core_solver_rows.size(); }
-
-    void fill_costs_and_x_for_first_stage_solver(unsigned original_number_of_columns);
-
-    void init_buffer(unsigned k, vector<T> & r);
-
-    void refactor();
-
-    void set_scaled_costs();
-public:
-    lp_primal_simplex(): m_core_solver(nullptr) {}
-
-    column_info<T> * get_or_create_column_info(unsigned column);
-
-    void set_status(lp_status status) {
-        this->m_status = status;
-    }
-
-    lp_status get_status() {
-        return this->m_status;
-    }
-
-    void fill_acceptable_values_for_x();
-
-
-    void set_zero_bound(bool * bound_is_set, T * bounds,  unsigned i);
-
-    void fill_costs_and_x_for_first_stage_solver_for_row(
-                                                         int row,
-                                                         unsigned & slack_var,
-                                                         unsigned & artificial);
-
-
-
-    
-    void set_core_solver_bounds();
-
-    void find_maximal_solution() override;
-
-    void fill_A_x_and_basis_for_stage_one_total_inf();
-
-    void fill_A_x_and_basis_for_stage_one_total_inf_for_row(unsigned row);
-
-    void solve_with_total_inf();
-
-
-    ~lp_primal_simplex() override;
-
-    bool bounds_hold(std::unordered_map<std::string, T> const & solution);
-
-    T get_row_value(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream * out);
-
-    bool row_constraint_holds(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream * out);
-
-    bool row_constraints_hold(std::unordered_map<std::string, T> const & solution);
-
-
-    T * get_array_from_map(std::unordered_map<std::string, T> const & solution);
-
-    bool solution_is_feasible(std::unordered_map<std::string, T> const & solution) {
-        return bounds_hold(solution) && row_constraints_hold(solution);
-    }
-
-    T get_column_value(unsigned column) const override {
-        return this->get_column_value_with_core_solver(column, m_core_solver);
-    }
-
-    T get_current_cost() const override;
-
-    
-};
-}
diff --git a/src/math/lp/lp_primal_simplex_def.h b/src/math/lp/lp_primal_simplex_def.h
deleted file mode 100644
index 7ffe819b2..000000000
--- a/src/math/lp/lp_primal_simplex_def.h
+++ /dev/null
@@ -1,367 +0,0 @@
-/*++
-Copyright (c) 2017 Microsoft Corporation
-
-Module Name:
-
-    <name>
-
-Abstract:
-
-    <abstract>
-
-Author:
-
-    Lev Nachmanson (levnach)
-
-Revision History:
-
-
---*/
-#pragma once
-
-#include <string>
-#include "util/vector.h"
-#include "math/lp/lp_primal_simplex.h"
-
-namespace lp {
-template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x_for_first_stage_solver(unsigned original_number_of_columns) {
-    unsigned slack_var = original_number_of_columns;
-    unsigned artificial = original_number_of_columns + this->m_slacks;
-
-    for (unsigned row = 0; row < this->row_count(); row++) {
-        fill_costs_and_x_for_first_stage_solver_for_row(row, slack_var, artificial);
-    }
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::init_buffer(unsigned k, vector<T> & r) {
-    for (unsigned i = 0; i < k; i++) {
-        r[i] = 0;
-    }
-    r[k] = 1;
-    for (unsigned i = this->row_count() -1; i > k; i--) {
-        r[i] = 0;
-    }
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::refactor() {
-    m_core_solver->init_lu();
-    if (m_core_solver->factorization()->get_status() != LU_status::OK) {
-        throw_exception("cannot refactor");
-    }
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::set_scaled_costs() {
-    unsigned j = this->number_of_core_structurals();
-    while (j-- > 0) {
-        this->set_scaled_cost(j);
-    }
-}
-
-template <typename T, typename X>     column_info<T> * lp_primal_simplex<T, X>::get_or_create_column_info(unsigned column) {
-    auto it = this->m_columns.find(column);
-    return (it == this->m_columns.end())? ( this->m_columns[column] = new column_info<T>) : it->second;
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::fill_acceptable_values_for_x() {
-    for (auto t : this->m_core_solver_columns_to_external_columns) {
-        this->m_x[t.first] = numeric_traits<T>::zero();
-    }
-}
-
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::set_zero_bound(bool * bound_is_set, T * bounds,  unsigned i) {
-    bound_is_set[i] = true;
-    bounds[i] = numeric_traits<T>::zero();
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::fill_costs_and_x_for_first_stage_solver_for_row(
-                                                                                                                int row,
-                                                                                                                unsigned & slack_var,
-                                                                                                                unsigned & artificial) {
-    lp_assert(row >= 0 && row < this->row_count());
-    auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
-    // we need to bring the program to the form Ax = b
-    T rs = this->m_b[row];
-    T artificial_cost =  - numeric_traits<T>::one();
-    switch (constraint.m_relation) {
-    case Equal: // no slack variable here
-        this->m_column_types[artificial] = column_type::lower_bound;
-        this->m_costs[artificial] = artificial_cost; // we are maximizing, so the artificial, which is non-negatiive, will be pushed to zero
-        this->m_basis[row] = artificial;
-        if (rs >= 0) {
-            (*this->m_A)(row, artificial) = numeric_traits<T>::one();
-            this->m_x[artificial] = rs;
-        } else {
-            (*this->m_A)(row, artificial) = - numeric_traits<T>::one();
-            this->m_x[artificial] = - rs;
-        }
-        artificial++;
-        break;
-
-    case Greater_or_equal:
-        this->m_column_types[slack_var] = column_type::lower_bound;
-        (*this->m_A)(row, slack_var) = - numeric_traits<T>::one();
-
-        if (rs > 0) {
-            lp_assert(numeric_traits<T>::is_zero(this->m_x[slack_var]));
-            // adding one artificial
-            this->m_column_types[artificial] = column_type::lower_bound;
-            (*this->m_A)(row, artificial) = numeric_traits<T>::one();
-            this->m_costs[artificial] = artificial_cost;
-            this->m_basis[row] = artificial;
-            this->m_x[artificial] = rs;
-            artificial++;
-        } else {
-            // we can put a slack_var into the basis, and atemplate <typename T, typename X> void lp_primal_simplex<T, X>::adding an artificial variable
-            this->m_basis[row] = slack_var;
-            this->m_x[slack_var] = - rs;
-        }
-        slack_var++;
-        break;
-    case Less_or_equal:
-        // introduce a non-negative slack variable
-        this->m_column_types[slack_var] = column_type::lower_bound;
-        (*this->m_A)(row, slack_var) = numeric_traits<T>::one();
-
-        if (rs < 0) {
-            // adding one artificial
-            lp_assert(numeric_traits<T>::is_zero(this->m_x[slack_var]));
-            this->m_column_types[artificial] = column_type::lower_bound;
-            (*this->m_A)(row, artificial) = - numeric_traits<T>::one();
-            this->m_costs[artificial] = artificial_cost;
-            this->m_x[artificial] = - rs;
-            this->m_basis[row] = artificial++;
-        } else {
-            // we can put slack_var into the basis, and atemplate <typename T, typename X> void lp_primal_simplex<T, X>::adding an artificial variable
-            this->m_basis[row] = slack_var;
-            this->m_x[slack_var] = rs;
-        }
-        slack_var++;
-        break;
-    }
-}
-
-
-
-
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::set_core_solver_bounds() {
-    unsigned total_vars = this->m_A->column_count() + this->m_slacks + this->m_artificials;
-    this->m_column_types.resize(total_vars);
-    this->m_upper_bounds.resize(total_vars);
-    for (auto cit : this->m_map_from_var_index_to_column_info) {
-        column_info<T> * ci = cit.second;
-        unsigned j = ci->get_column_index();
-        if (!is_valid(j))
-            continue; // the variable is not mapped to a column
-        switch (this->m_column_types[j] = ci->get_column_type()){
-        case column_type::fixed:
-            this->m_upper_bounds[j] = numeric_traits<T>::zero();
-            break;
-        case column_type::boxed:
-            this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j];
-            break;
-
-        default: break; // do nothing
-        }
-    }
-}
-
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::find_maximal_solution() {
-    if (this->problem_is_empty()) {
-        this->m_status = lp_status::EMPTY;
-        return;
-    }
-
-    this->cleanup();
-    this->fill_matrix_A_and_init_right_side();
-    if (this->m_status == lp_status::INFEASIBLE) {
-        return;
-    }
-    this->m_x.resize(this->m_A->column_count());
-    this->fill_m_b();
-    this->scale();
-    fill_acceptable_values_for_x();
-    this->count_slacks_and_artificials();
-    set_core_solver_bounds();
-    solve_with_total_inf();
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::fill_A_x_and_basis_for_stage_one_total_inf() {
-    for (unsigned row = 0; row < this->row_count(); row++)
-        fill_A_x_and_basis_for_stage_one_total_inf_for_row(row);
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::fill_A_x_and_basis_for_stage_one_total_inf_for_row(unsigned row) {
-    lp_assert(row < this->row_count());
-    auto ext_row_it = this->m_core_solver_rows_to_external_rows.find(row);
-    lp_assert(ext_row_it != this->m_core_solver_rows_to_external_rows.end());
-    unsigned ext_row = ext_row_it->second;
-    auto constr_it = this->m_constraints.find(ext_row);
-    lp_assert(constr_it != this->m_constraints.end());
-    auto & constraint = constr_it->second;
-    unsigned j = this->m_A->column_count(); // j is a slack variable
-    this->m_A->add_column();
-    // we need to bring the program to the form Ax = b
-    this->m_basis[row] = j;
-    switch (constraint.m_relation) {
-    case Equal:
-        this->m_x[j] = this->m_b[row];
-        (*this->m_A)(row, j) = numeric_traits<T>::one();
-        this->m_column_types[j] = column_type::fixed;
-        this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type<X>();
-        break;
-
-    case Greater_or_equal:
-        this->m_x[j] = - this->m_b[row];
-        (*this->m_A)(row, j) = - numeric_traits<T>::one();
-        this->m_column_types[j] = column_type::lower_bound;
-        this->m_upper_bounds[j] = zero_of_type<X>();
-        break;
-    case Less_or_equal:
-        this->m_x[j] = this->m_b[row];
-        (*this->m_A)(row, j) = numeric_traits<T>::one();
-        this->m_column_types[j] = column_type::lower_bound;
-        this->m_upper_bounds[j] = m_lower_bounds[j] = zero_of_type<X>();
-        break;
-    default:
-        lp_unreachable();
-    }
-}
-
-template <typename T, typename X> void lp_primal_simplex<T, X>::solve_with_total_inf() {
-    int total_vars = this->m_A->column_count() + this->row_count();
-    if (total_vars == 0) {
-        this->m_status = lp_status::OPTIMAL;
-        return;
-    }
-    m_lower_bounds.clear();
-    m_lower_bounds.resize(total_vars, zero_of_type<X>());  // low bounds are shifted ot zero
-    this->m_x.resize(total_vars, numeric_traits<T>::zero());
-    this->m_basis.resize(this->row_count());
-    this->m_costs.clear();
-    this->m_costs.resize(total_vars, zero_of_type<T>());
-    fill_A_x_and_basis_for_stage_one_total_inf();
-    if (this->m_settings.get_message_ostream() != nullptr)
-        this->print_statistics_on_A(*this->m_settings.get_message_ostream());
-    set_scaled_costs();
-
-    m_core_solver = new lp_primal_core_solver<T, X>(*this->m_A,
-                                                    this->m_b,
-                                                    this->m_x,
-                                                    this->m_basis,
-                                                    this->m_nbasis,
-                                                    this->m_heading,
-                                                    this->m_costs,
-                                                    this->m_column_types,
-                                                    m_lower_bounds,
-                                                    this->m_upper_bounds,
-                                                    this->m_settings, *this);
-    m_core_solver->solve();
-    this->set_status(m_core_solver->get_status());
-    this->m_total_iterations = m_core_solver->total_iterations();   
-}
-
-
-template <typename T, typename X> lp_primal_simplex<T, X>::~lp_primal_simplex() {
-    delete m_core_solver;
-}
-
-template <typename T, typename X> bool lp_primal_simplex<T, X>::bounds_hold(std::unordered_map<std::string, T> const & solution) {
-    for (auto it : this->m_map_from_var_index_to_column_info) {
-        auto sol_it = solution.find(it.second->get_name());
-        if (sol_it == solution.end()) {
-            std::stringstream s;
-            s << "cannot find column " << it.first << " in solution";
-            throw_exception(s.str() );
-        }
-
-        if (!it.second->bounds_hold(sol_it->second)) {
-            it.second->bounds_hold(sol_it->second);
-            return false;
-        }
-    }
-    return true;
-}
-
-template <typename T, typename X> T lp_primal_simplex<T, X>::get_row_value(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream * out) {
-    auto it = this->m_A_values.find(i);
-    if (it == this->m_A_values.end()) {
-        std::stringstream s;
-        s << "cannot find row " << i;
-        throw_exception(s.str() );
-    }
-    T ret = numeric_traits<T>::zero();
-    for (auto & pair : it->second) {
-        auto cit = this->m_map_from_var_index_to_column_info.find(pair.first);
-        lp_assert(cit != this->m_map_from_var_index_to_column_info.end());
-        column_info<T> * ci = cit->second;
-        auto sol_it = solution.find(ci->get_name());
-        lp_assert(sol_it != solution.end());
-        T column_val = sol_it->second;
-        if (out != nullptr) {
-            (*out) << pair.second << "(" << ci->get_name() << "=" << column_val << ") ";
-        }
-        ret += pair.second * column_val;
-    }
-    if (out != nullptr) {
-        (*out) << " = " << ret << std::endl;
-    }
-    return ret;
-}
-
-template <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraint_holds(unsigned i, std::unordered_map<std::string, T> const & solution, std::ostream *out) {
-    T row_val = get_row_value(i, solution, out);
-    auto & constraint = this->m_constraints[i];
-    T rs = constraint.m_rs;
-    bool print = out != nullptr;
-    switch (constraint.m_relation) {
-    case Equal:
-        if (fabs(numeric_traits<T>::get_double(row_val - rs)) > 0.00001) {
-            if (print) {
-                (*out) << "should be = " << rs << std::endl;
-            }
-            return false;
-        }
-        return true;
-    case Greater_or_equal:
-        if (numeric_traits<T>::get_double(row_val - rs) < -0.00001) {
-            if (print) {
-                (*out) << "should be >= " << rs << std::endl;
-            }
-            return false;
-        }
-        return true;;
-
-    case Less_or_equal:
-        if (numeric_traits<T>::get_double(row_val - rs) > 0.00001) {
-            if (print) {
-                (*out) << "should be <= " << rs << std::endl;
-            }
-            return false;
-        }
-        return true;;
-    }
-    lp_unreachable();
-    return false; // it is unreachable
-}
-
-template <typename T, typename X> bool lp_primal_simplex<T, X>::row_constraints_hold(std::unordered_map<std::string, T> const & solution) {
-    for (auto it : this->m_A_values) {
-        if (!row_constraint_holds(it.first, solution, nullptr)) {
-            row_constraint_holds(it.first, solution, nullptr);
-            return false;
-        }
-    }
-    return true;
-}
-
-template <typename T, typename X> T lp_primal_simplex<T, X>::get_current_cost() const {
-    T ret = numeric_traits<T>::zero();
-    for (auto it : this->m_map_from_var_index_to_column_info) {
-        ret += this->get_column_cost_value(it.first, it.second);
-    }
-    return ret;
-}
-}
diff --git a/src/sat/smt/arith_sls.h b/src/sat/smt/arith_sls.h
index f50441871..b0b4fb48b 100644
--- a/src/sat/smt/arith_sls.h
+++ b/src/sat/smt/arith_sls.h
@@ -20,7 +20,6 @@ Author:
 #include "ast/ast_trail.h"
 #include "ast/arith_decl_plugin.h"
 #include "math/lp/lp_solver.h"
-#include "math/lp/lp_primal_simplex.h"
 #include "math/lp/indexed_value.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/nla_solver.h"
diff --git a/src/sat/smt/arith_solver.h b/src/sat/smt/arith_solver.h
index eac73c719..3152485ec 100644
--- a/src/sat/smt/arith_solver.h
+++ b/src/sat/smt/arith_solver.h
@@ -20,7 +20,6 @@ Author:
 #include "ast/ast_trail.h"
 #include "ast/arith_decl_plugin.h"
 #include "math/lp/lp_solver.h"
-#include "math/lp/lp_primal_simplex.h"
 #include "math/lp/indexed_value.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/nla_solver.h"
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 4180f9ec3..a654366a2 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -20,7 +20,6 @@
   --*/
 #include "util/stopwatch.h"
 #include "math/lp/lp_solver.h"
-#include "math/lp/lp_primal_simplex.h"
 #include "math/lp/indexed_value.h"
 #include "math/lp/lar_solver.h"
 #include "math/lp/nla_solver.h"
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index c2b297e6b..78abf1a6f 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -33,7 +33,6 @@
 #include <stdlib.h>
 #include <utility>
 #include "math/lp/lp_utils.h"
-#include "math/lp/lp_primal_simplex.h"
 #include "test/lp/smt_reader.h"
 #include "math/lp/binary_heap_priority_queue.h"
 #include "test/lp/argument_parser.h"
diff --git a/src/test/lp/smt_reader.h b/src/test/lp/smt_reader.h
index 6060370a8..75edb23b9 100644
--- a/src/test/lp/smt_reader.h
+++ b/src/test/lp/smt_reader.h
@@ -23,7 +23,6 @@ Revision History:
 #include <string>
 #include <vector>
 #include <unordered_map>
-#include "math/lp/lp_primal_simplex.h"
 #include "math/lp/lar_solver.h"
 #include <iostream>
 #include <fstream>