diff --git a/README-CMake.md b/README-CMake.md
index 9289d40f1..abd2d3c87 100644
--- a/README-CMake.md
+++ b/README-CMake.md
@@ -1,4 +1,4 @@
-# Z3's CMake build system
+    # Z3's CMake build system
 
 [CMake](https://cmake.org/) is a "meta build system" that reads a description
 of the project written in the ``CMakeLists.txt`` files and emits a build
diff --git a/src/shell/lp_frontend.cpp b/src/shell/lp_frontend.cpp
index 9fad426a0..8bed24c83 100644
--- a/src/shell/lp_frontend.cpp
+++ b/src/shell/lp_frontend.cpp
@@ -80,7 +80,6 @@ void run_solver(lp_params & params, char const * mps_file_name) {
     solver->settings().set_message_ostream(&std::cout);
     solver->settings().report_frequency = params.rep_freq();
     solver->settings().print_statistics = params.print_stats();
-    solver->settings().presolve_with_double_solver_for_lar = params.presolve_with_dbl();
     solver->find_maximal_solution();
 
     *(solver->settings().get_message_ostream()) << "status is " << lp_status_to_string(solver->get_status()) << std::endl;
diff --git a/src/smt/theory_lra.cpp b/src/smt/theory_lra.cpp
index 4a602f808..4cff79cf2 100644
--- a/src/smt/theory_lra.cpp
+++ b/src/smt/theory_lra.cpp
@@ -148,8 +148,6 @@ namespace smt {
         arith_eq_adapter     m_arith_eq_adapter;
 
         vector<rational>    m_columns;
-        int m_print_counter = 0;
-
         // temporary values kept during internalization
         struct internalize_state {
             expr_ref_vector     m_terms;                     
@@ -287,7 +285,6 @@ namespace smt {
             m_theory_var2var_index.reset();
             m_solver->settings().set_resource_limit(m_resource_limit);
             m_solver->settings().simplex_strategy() = static_cast<lean::simplex_strategy_enum>(lp.simplex_strategy());
-            m_solver->settings().presolve_with_double_solver_for_lar = lp.presolve_with_dbl();
             reset_variable_values();
             m_solver->settings().bound_propagation() = BP_NONE != propagation_mode();
             m_solver->set_propagate_bounds_on_pivoted_rows_mode(lp.bprop_on_pivoted_rows());
@@ -1985,7 +1982,6 @@ namespace smt {
         typedef pair_hash<obj_hash<rational>, bool_hash> value_sort_pair_hash;
         typedef map<value_sort_pair, theory_var, value_sort_pair_hash, default_eq<value_sort_pair> > value2var;
         value2var               m_fixed_var_table;
-        const lean::constraint_index null_index = static_cast<lean::constraint_index>(-1);
 
         void propagate_eqs(lean::var_index vi, lean::constraint_index ci, lean::lconstraint_kind k, lp::bound& b) {
             if (propagate_eqs()) {
@@ -2029,10 +2025,10 @@ namespace smt {
             lean::var_index ti = m_solver->adjust_term_index(vi);
             auto& vec = is_lower ? m_lower_terms : m_upper_terms;
             if (vec.size() <= ti) {
-                vec.resize(ti + 1, constraint_bound(null_index, rational()));
+                vec.resize(ti + 1, constraint_bound(UINT_MAX, rational()));
             }
             constraint_bound& b = vec[ti];
-            if (b.first == null_index || (is_lower? b.second < v : b.second > v)) {
+            if (b.first == UINT_MAX || (is_lower? b.second < v : b.second > v)) {
                 ctx().push_trail(vector_value_trail<context, constraint_bound>(vec, ti));
                 b.first = ci;
                 b.second = v;
@@ -2052,7 +2048,7 @@ namespace smt {
                 rational val;
                 TRACE("arith", tout << vi << " " << v << "\n";);
                 if (v != null_theory_var && a.is_numeral(get_owner(v), val) && bound == val) {
-                    ci = null_constraint_index;
+                    ci = UINT_MAX;
                     return bound == val;
                 }
 
@@ -2060,7 +2056,7 @@ namespace smt {
                 if (vec.size() > ti) {
                     constraint_bound& b = vec[ti];
                     ci = b.first;
-                    return ci != null_index && bound == b.second;
+                    return ci != UINT_MAX && bound == b.second;
                 }
                 else {
                     return false;
@@ -2144,22 +2140,10 @@ namespace smt {
             if (m_solver->A_r().row_count() > m_stats.m_max_rows)
                 m_stats.m_max_rows = m_solver->A_r().row_count();
             TRACE("arith_verbose", display(tout););
-            bool print = false && m_print_counter++ % 1000 == 0;
-            stopwatch sw;
-            if (print) {
-                sw.start();
-            }
             lean::lp_status status = m_solver->find_feasible_solution();
-            if (print) {
-                sw.stop();
-            }
             m_stats.m_num_iterations = m_solver->settings().st().m_total_iterations;
             m_stats.m_num_factorizations = m_solver->settings().st().m_num_factorizations;
             m_stats.m_need_to_solve_inf = m_solver->settings().st().m_need_to_solve_inf;
-            if (print) {
-                IF_VERBOSE(0, verbose_stream() << status << " " << sw.get_seconds() << " " << m_stats.m_num_iterations << " " << m_print_counter << "\n";);
-            }
-            //m_stats.m_num_iterations_with_no_progress += m_solver->settings().st().m_iters_with_no_cost_growing;
 
             switch (status) {
             case lean::lp_status::INFEASIBLE:
@@ -2182,10 +2166,11 @@ namespace smt {
         literal_vector      m_core;
         svector<enode_pair> m_eqs;
         vector<parameter>   m_params;
-        lean::constraint_index const null_constraint_index = UINT_MAX;
+
+        // lean::constraint_index const null_constraint_index = UINT_MAX; // not sure what a correct fix is
 
         void set_evidence(lean::constraint_index idx) {
-            if (idx == null_constraint_index) {
+            if (idx == UINT_MAX) {
                 return;
             }
             switch (m_constraint_sources[idx]) {
diff --git a/src/test/lp.cpp b/src/test/lp.cpp
index 7829d81f5..f1df5b05e 100644
--- a/src/test/lp.cpp
+++ b/src/test/lp.cpp
@@ -1077,7 +1077,7 @@ bool get_double_from_args_parser(const char * option, argument_parser & args_par
 
 void update_settings(argument_parser & args_parser, lp_settings& settings) {
     unsigned n;
-    settings.m_simplex_strategy = simplex_strategy_enum::no_tableau;
+    settings.m_simplex_strategy = simplex_strategy_enum::lu;
     if (get_int_from_args_parser("--rep_frq", args_parser, n))
         settings.report_frequency = n;
     else
diff --git a/src/test/smt_reader.h b/src/test/smt_reader.h
index 481985ed4..714e734ad 100644
--- a/src/test/smt_reader.h
+++ b/src/test/smt_reader.h
@@ -56,10 +56,11 @@ namespace lean {
         struct formula_constraint {
             lconstraint_kind m_kind;
             std::vector<std::pair<mpq, std::string>> m_coeffs;
-            mpq m_right_side = numeric_traits<mpq>::zero();
+            mpq m_right_side;
             void add_pair(mpq c, std::string name) {
                 m_coeffs.push_back(make_pair(c, name));
             }
+            formula_constraint() : m_right_side(numeric_traits<mpq>::zero()) {}
         };
 
         lisp_elem m_formula_lisp_elem;
@@ -69,9 +70,11 @@ namespace lean {
         std::string m_file_name;
         std::ifstream m_file_stream;
         std::string m_line;
-        bool m_is_OK = true;
-        unsigned m_line_number = 0;
-        smt_reader(std::string file_name): 
+        bool m_is_OK;
+        unsigned m_line_number;
+        smt_reader(std::string file_name):
+            m_is_OK(true),
+            m_line_number(0),
             m_file_name(file_name), m_file_stream(file_name) {
         }
 
diff --git a/src/util/lp/binary_heap_priority_queue.h b/src/util/lp/binary_heap_priority_queue.h
index 1a1f89592..a6206948c 100644
--- a/src/util/lp/binary_heap_priority_queue.h
+++ b/src/util/lp/binary_heap_priority_queue.h
@@ -16,8 +16,7 @@ class binary_heap_priority_queue {
     // indexing for A starts from 1
     vector<unsigned> m_heap; // keeps the elements of the queue
     vector<int> m_heap_inverse; // o = m_heap[m_heap_inverse[o]]
-    unsigned m_heap_size = 0;
-
+    unsigned m_heap_size;
     // is is the child place in heap
     void swap_with_parent(unsigned i);
     void put_at(unsigned i, unsigned h);
@@ -29,7 +28,7 @@ public:
 public:
     void remove(unsigned o);
     unsigned size() const { return m_heap_size; }
-    binary_heap_priority_queue(): m_heap(1) {} // the empty constructror
+    binary_heap_priority_queue(): m_heap(1), m_heap_size(0) {} // the empty constructror
     // n is the initial queue capacity.
     // The capacity will be enlarged two times automatically if needed
     binary_heap_priority_queue(unsigned n);
diff --git a/src/util/lp/binary_heap_priority_queue.hpp b/src/util/lp/binary_heap_priority_queue.hpp
index 2ad65c536..440b45b02 100644
--- a/src/util/lp/binary_heap_priority_queue.hpp
+++ b/src/util/lp/binary_heap_priority_queue.hpp
@@ -83,7 +83,8 @@ template <typename T> void binary_heap_priority_queue<T>::remove(unsigned o) {
 template <typename T> binary_heap_priority_queue<T>::binary_heap_priority_queue(unsigned n) :
     m_priorities(n),
     m_heap(n + 1), // because the indexing for A starts from 1
-    m_heap_inverse(n, -1)
+    m_heap_inverse(n, -1),
+    m_heap_size(0)
 { }
 
 
diff --git a/src/util/lp/binary_heap_upair_queue.h b/src/util/lp/binary_heap_upair_queue.h
index d5df7affc..26cfd5532 100644
--- a/src/util/lp/binary_heap_upair_queue.h
+++ b/src/util/lp/binary_heap_upair_queue.h
@@ -20,8 +20,8 @@ template <typename  T>
 class binary_heap_upair_queue {
     binary_heap_priority_queue<T> m_q;
     std::unordered_map<upair, unsigned> m_pairs_to_index;
-    vector<upair> m_pairs; // inverse to index
-    vector<unsigned> m_available_spots;
+    svector<upair> m_pairs; // inverse to index
+    svector<unsigned> m_available_spots;
 public:
     binary_heap_upair_queue(unsigned size);
 
diff --git a/src/util/lp/bound_analyzer_on_row.h b/src/util/lp/bound_analyzer_on_row.h
index 7987c89e9..bc0b72eee 100644
--- a/src/util/lp/bound_analyzer_on_row.h
+++ b/src/util/lp/bound_analyzer_on_row.h
@@ -19,9 +19,9 @@ class bound_analyzer_on_row {
     
     linear_combination_iterator<mpq> & m_it;
     unsigned m_row_or_term_index;
-    int m_column_of_u = -1; // index of an unlimited from above monoid
+    int m_column_of_u; // index of an unlimited from above monoid
     // -1 means that such a value is not found, -2 means that at least two of such monoids were found
-    int m_column_of_l = -1; // index of an unlimited from below monoid
+    int m_column_of_l; // index of an unlimited from below monoid
     impq m_rs;
     bound_propagator & m_bp;
 public :
@@ -36,7 +36,9 @@ public :
         m_it(it),
         m_row_or_term_index(row_or_term_index),
         m_rs(rs),
-        m_bp(bp)
+        m_bp(bp),
+        m_column_of_u(-1),
+        m_column_of_l(-1)
     {}
 
 
diff --git a/src/util/lp/column_info.h b/src/util/lp/column_info.h
index 7e44ccf43..ed916757d 100644
--- a/src/util/lp/column_info.h
+++ b/src/util/lp/column_info.h
@@ -15,16 +15,16 @@ inline bool is_valid(unsigned j) { return static_cast<int>(j) >= 0;}
 template <typename T>
 class column_info {
     std::string m_name;
-    bool m_low_bound_is_set = false;
-    bool m_low_bound_is_strict = false;
-    bool m_upper_bound_is_set = false;
-    bool m_upper_bound_is_strict = false;
+    bool m_low_bound_is_set;
+    bool m_low_bound_is_strict;
+    bool m_upper_bound_is_set;
+    bool m_upper_bound_is_strict;
     T m_low_bound;
     T m_upper_bound;
-    T m_cost = numeric_traits<T>::zero();
+    T m_cost;
     T m_fixed_value;
-    bool m_is_fixed = false;
-    unsigned m_column_index = static_cast<unsigned>(-1);
+    bool m_is_fixed;
+    unsigned m_column_index;
 public:
     bool operator==(const column_info & c) const {
         return     m_name == c.m_name &&
@@ -44,9 +44,24 @@ public:
         m_column_index = j;
     }
     // the default constructor
-    column_info() {}
-
-    column_info(unsigned column_index) : m_column_index(column_index) {
+    column_info():
+        m_low_bound_is_set(false),
+        m_low_bound_is_strict(false),
+        m_upper_bound_is_set (false),
+        m_upper_bound_is_strict (false),
+        m_is_fixed(false),
+        m_cost(numeric_traits<T>::zero()),
+        m_column_index(static_cast<unsigned>(-1))
+    {}
+    
+    column_info(unsigned column_index) :
+        m_low_bound_is_set(false),
+        m_low_bound_is_strict(false),
+        m_upper_bound_is_set (false),
+        m_upper_bound_is_strict (false),
+        m_is_fixed(false),
+        m_cost(numeric_traits<T>::zero()),
+        m_column_index(column_index) {
     }
 
     column_info(const column_info & ci) {
diff --git a/src/util/lp/conversion_helper.h b/src/util/lp/conversion_helper.h
new file mode 100644
index 000000000..ea88a98bf
--- /dev/null
+++ b/src/util/lp/conversion_helper.h
@@ -0,0 +1,44 @@
+/*
+  Copyright (c) 2013 Microsoft Corporation. All rights reserved.
+  Released under Apache 2.0 license as described in the file LICENSE.
+
+  Author: Lev Nachmanson
+*/
+#pragma once
+namespace lean {
+template <typename V>
+struct conversion_helper {
+    static V get_low_bound(const column_info<mpq> & ci) {
+        return V(ci.get_low_bound(), ci.low_bound_is_strict()? 1 : 0);
+    }
+
+    static V get_upper_bound(const column_info<mpq> & ci) {
+        return V(ci.get_upper_bound(), ci.upper_bound_is_strict()? -1 : 0);
+    }
+};
+
+template<>
+struct conversion_helper <double> {
+    static double conversion_helper <double>::get_upper_bound(const column_info<mpq> & ci) {
+        if (!ci.upper_bound_is_strict())
+            return ci.get_upper_bound().get_double();
+        double eps = 0.00001;
+        if (!ci.low_bound_is_set())
+            return ci.get_upper_bound().get_double() - eps;
+        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
+        return ci.get_upper_bound().get_double() - eps;
+    }
+
+    static double get_low_bound(const column_info<mpq> & ci) {
+        if (!ci.low_bound_is_strict())
+            return ci.get_low_bound().get_double();
+        double eps = 0.00001;
+        if (!ci.upper_bound_is_set())
+            return ci.get_low_bound().get_double() + eps;
+        eps = std::min((ci.get_upper_bound() - ci.get_low_bound()).get_double() / 1000, eps);
+        return ci.get_low_bound().get_double() + eps;
+    }
+
+};
+
+}
diff --git a/src/util/lp/core_solver_pretty_printer.h b/src/util/lp/core_solver_pretty_printer.h
index 6dbeb28cf..2a3a14b31 100644
--- a/src/util/lp/core_solver_pretty_printer.h
+++ b/src/util/lp/core_solver_pretty_printer.h
@@ -16,7 +16,6 @@ template <typename T, typename X> class lp_core_solver_base; // forward definiti
 template <typename T, typename X>
 class core_solver_pretty_printer {
     std::ostream & m_out;
-    template<typename A> using vector = vector<A>;
     typedef std::string string;
     lp_core_solver_base<T, X> & m_core_solver;
     vector<unsigned> m_column_widths;
@@ -34,15 +33,15 @@ class core_solver_pretty_printer {
     std::string m_cost_title;
     std::string m_basis_heading_title;
     std::string m_x_title;
-    std::string m_low_bounds_title = "low";
-    std::string m_upp_bounds_title = "upp";
-    std::string m_exact_norm_title = "exact cn";
-    std::string m_approx_norm_title = "approx cn";
+    std::string m_low_bounds_title;
+    std::string m_upp_bounds_title;
+    std::string m_exact_norm_title;
+    std::string m_approx_norm_title;
 
 
     unsigned ncols() { return m_core_solver.m_A.column_count(); }
     unsigned nrows() { return m_core_solver.m_A.row_count(); }
-    unsigned m_artificial_start = std::numeric_limits<unsigned>::max();
+    unsigned m_artificial_start;
     indexed_vector<T> m_w_buff;
     indexed_vector<T> m_ed_buff;
     vector<T> m_exact_column_norms;
diff --git a/src/util/lp/core_solver_pretty_printer.hpp b/src/util/lp/core_solver_pretty_printer.hpp
index aa963b220..786b8b3a1 100644
--- a/src/util/lp/core_solver_pretty_printer.hpp
+++ b/src/util/lp/core_solver_pretty_printer.hpp
@@ -23,6 +23,12 @@ core_solver_pretty_printer<T, X>::core_solver_pretty_printer(lp_core_solver_base
     m_rs(ncols(), zero_of_type<X>()),
     m_w_buff(core_solver.m_w),
     m_ed_buff(core_solver.m_ed) {
+    m_low_bounds_title = "low";
+    m_upp_bounds_title = "upp";
+    m_exact_norm_title = "exact cn";
+    m_approx_norm_title = "approx cn";
+    m_artificial_start = std::numeric_limits<unsigned>::max();
+
     m_column_widths.resize(core_solver.m_A.column_count(), 0),
     init_m_A_and_signs();
     init_costs();
diff --git a/src/util/lp/init_lar_solver.h b/src/util/lp/init_lar_solver.h
new file mode 100644
index 000000000..3c74e622d
--- /dev/null
+++ b/src/util/lp/init_lar_solver.h
@@ -0,0 +1,573 @@
+/*
+  Copyright (c) 2017 Microsoft Corporation
+  Author: Lev Nachmanson
+*/
+
+// here we are inside lean::lar_solver class
+
+bool strategy_is_undecided() const {
+    return m_settings.simplex_strategy() == simplex_strategy_enum::undecided;
+}
+
+var_index add_var(unsigned ext_j) {
+    var_index i;
+    lean_assert (ext_j < m_terms_start_index); 
+
+    if (ext_j >= m_terms_start_index)
+        throw 0; // todo : what is the right way to exit?
+            
+    if (try_get_val(m_ext_vars_to_columns, ext_j, i)) {
+        return i;
+    }
+    lean_assert(m_vars_to_ul_pairs.size() == A_r().column_count());
+    i = A_r().column_count();
+    m_vars_to_ul_pairs.push_back (ul_pair(static_cast<unsigned>(-1)));
+    add_non_basic_var_to_core_fields(ext_j);
+    lean_assert(sizes_are_correct());
+    return i;
+}
+
+void register_new_ext_var_index(unsigned ext_v) {
+    lean_assert(!contains(m_ext_vars_to_columns, ext_v));
+    unsigned j = m_ext_vars_to_columns.size();
+    m_ext_vars_to_columns[ext_v] = j;
+    lean_assert(m_columns_to_ext_vars_or_term_indices.size() == j);
+    m_columns_to_ext_vars_or_term_indices.push_back(ext_v);
+}
+
+void add_non_basic_var_to_core_fields(unsigned ext_j) {
+    register_new_ext_var_index(ext_j);
+    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
+    m_columns_with_changed_bound.increase_size_by_one();
+    add_new_var_to_core_fields_for_mpq(false);
+    if (use_lu())
+        add_new_var_to_core_fields_for_doubles(false);
+}
+
+void add_new_var_to_core_fields_for_doubles(bool register_in_basis) {
+    unsigned j = A_d().column_count();
+    A_d().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
+    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
+    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_d().add_row();
+        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
+        m_mpq_lar_core_solver.m_d_basis.push_back(j);
+    }else {
+        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
+    }
+}
+
+void add_new_var_to_core_fields_for_mpq(bool register_in_basis) {
+    unsigned j = A_r().column_count();
+    A_r().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_r_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_r_low_bounds.size() == j && m_mpq_lar_core_solver.m_r_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_r_x.resize(j + 1);
+    m_mpq_lar_core_solver.m_r_low_bounds.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_upper_bounds.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_solver.m_inf_set.increase_size_by_one();
+    m_mpq_lar_core_solver.m_r_solver.m_costs.resize(j + 1);
+    m_mpq_lar_core_solver.m_r_solver.m_d.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_r_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_r().add_row();
+        m_mpq_lar_core_solver.m_r_heading.push_back(m_mpq_lar_core_solver.m_r_basis.size());
+        m_mpq_lar_core_solver.m_r_basis.push_back(j);
+        if (m_settings.bound_propagation())
+            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
+    } else {
+        m_mpq_lar_core_solver.m_r_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_r_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_r_nbasis.push_back(j);
+    }
+}
+
+
+var_index add_term_undecided(const vector<std::pair<mpq, var_index>> & coeffs,
+                             const mpq &m_v) {
+    m_terms.push_back(new lar_term(coeffs, m_v));
+    m_orig_terms.push_back(new lar_term(coeffs, m_v));
+    return m_terms_start_index +  m_terms.size() - 1;
+}
+
+// terms
+var_index add_term(const vector<std::pair<mpq, var_index>> & coeffs,
+                   const mpq &m_v) {
+    if (strategy_is_undecided())
+        return add_term_undecided(coeffs, m_v);
+
+    m_terms.push_back(new lar_term(coeffs, m_v));
+    m_orig_terms.push_back(new lar_term(coeffs, m_v));
+    unsigned adjusted_term_index = m_terms.size() - 1;
+	var_index ret = m_terms_start_index + adjusted_term_index;
+    if (use_tableau() && !coeffs.empty()) {
+        add_row_for_term(m_orig_terms.back(), ret);
+        if (m_settings.bound_propagation())
+            m_rows_with_changed_bounds.insert(A_r().row_count() - 1);
+    } 
+    lean_assert(m_ext_vars_to_columns.size() == A_r().column_count());
+    return ret;
+}
+
+void add_row_for_term(const lar_term * term, unsigned term_ext_index) {
+    lean_assert(sizes_are_correct());
+    add_row_from_term_no_constraint(term, term_ext_index);
+    lean_assert(sizes_are_correct());
+}
+
+void add_row_from_term_no_constraint(const lar_term * term, unsigned term_ext_index) {
+    register_new_ext_var_index(term_ext_index);
+    // j will be a new variable
+	unsigned j = A_r().column_count();
+    ul_pair ul(j);
+    m_vars_to_ul_pairs.push_back(ul);
+    add_basic_var_to_core_fields();
+    if (use_tableau()) {
+        auto it = iterator_on_term_with_basis_var(*term, j);
+        A_r().fill_last_row_with_pivoting(it,
+                                          m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
+        m_mpq_lar_core_solver.m_r_solver.m_b.resize(A_r().column_count(), zero_of_type<mpq>());
+    } else {
+        fill_last_row_of_A_r(A_r(), term);
+    }
+    m_mpq_lar_core_solver.m_r_x[j] = get_basic_var_value_from_row_directly(A_r().row_count() - 1);
+    if (use_lu())
+        fill_last_row_of_A_d(A_d(), term);
+}
+
+void add_basic_var_to_core_fields() {
+    bool use_lu = m_mpq_lar_core_solver.need_to_presolve_with_double_solver();
+    lean_assert(!use_lu || A_r().column_count() == A_d().column_count());
+    m_mpq_lar_core_solver.m_column_types.push_back(column_type::free_column);
+    m_columns_with_changed_bound.increase_size_by_one();
+    m_rows_with_changed_bounds.increase_size_by_one();
+    add_new_var_to_core_fields_for_mpq(true);
+    if (use_lu)
+        add_new_var_to_core_fields_for_doubles(true);
+}
+
+constraint_index add_var_bound(var_index j, lconstraint_kind kind, const mpq & right_side)  {
+	constraint_index ci = m_constraints.size();
+    if (!is_term(j)) { // j is a var
+        auto vc = new lar_var_constraint(j, kind, right_side);
+        m_constraints.push_back(vc);
+        update_column_type_and_bound(j, kind, right_side, ci);
+    } else {
+        add_var_bound_on_constraint_for_term(j, kind, right_side, ci);
+    }
+    lean_assert(sizes_are_correct());
+    return ci;
+}
+
+void update_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_index) {
+    switch(m_mpq_lar_core_solver.m_column_types[j]) {
+    case column_type::free_column:
+        update_free_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::boxed:
+        update_boxed_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::low_bound:
+        update_low_bound_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::upper_bound:
+        update_upper_bound_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    case column_type::fixed:
+        update_fixed_column_type_and_bound(j, kind, right_side, constr_index);
+        break;
+    default:
+        lean_assert(false); // cannot be here
+    }
+}
+
+void add_var_bound_on_constraint_for_term(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(is_term(j));
+    unsigned adjusted_term_index = adjust_term_index(j);
+    unsigned term_j;
+    if (try_get_val(m_ext_vars_to_columns, j, term_j)) {
+        mpq rs = right_side - m_orig_terms[adjusted_term_index]->m_v;
+        m_constraints.push_back(new lar_term_constraint(m_orig_terms[adjusted_term_index], kind, right_side));
+        update_column_type_and_bound(term_j, kind, rs, ci);
+    }
+    else {
+        add_constraint_from_term_and_create_new_column_row(j, m_orig_terms[adjusted_term_index], kind, right_side);
+    }
+}
+
+
+void add_constraint_from_term_and_create_new_column_row(unsigned term_j, const lar_term* term,
+                                                        lconstraint_kind kind, const mpq & right_side) {
+
+    add_row_from_term_no_constraint(term, term_j);
+    unsigned j = A_r().column_count() - 1;
+    update_column_type_and_bound(j, kind, right_side - term->m_v, m_constraints.size());
+    m_constraints.push_back(new lar_term_constraint(term, kind, right_side));
+    lean_assert(A_r().column_count() == m_mpq_lar_core_solver.m_r_solver.m_costs.size());
+}
+
+void decide_on_strategy_and_adjust_initial_state() {
+	lean_assert(strategy_is_undecided());
+	if (m_vars_to_ul_pairs.size() > m_settings.column_number_threshold_for_using_lu_in_lar_solver) {
+        m_settings.simplex_strategy() = simplex_strategy_enum::lu;
+    } else {
+        m_settings.simplex_strategy() = simplex_strategy_enum::tableau_rows; // todo: when to switch to tableau_costs?
+    }
+    adjust_initial_state();
+}
+
+void adjust_initial_state() {
+    switch (m_settings.simplex_strategy()) {
+    case simplex_strategy_enum::lu:
+        adjust_initial_state_for_lu();
+        break;
+    case simplex_strategy_enum::tableau_rows:
+        adjust_initial_state_for_tableau_rows();
+        break;
+    case simplex_strategy_enum::tableau_costs:
+        lean_assert(false); // not implemented
+    }
+}
+
+void adjust_initial_state_for_lu() {
+    copy_from_mpq_matrix(A_d());
+	unsigned n = A_d().column_count();
+	m_mpq_lar_core_solver.m_d_x.resize(n);
+	m_mpq_lar_core_solver.m_d_low_bounds.resize(n);
+	m_mpq_lar_core_solver.m_d_upper_bounds.resize(n);
+	m_mpq_lar_core_solver.m_d_heading = m_mpq_lar_core_solver.m_r_heading;
+	m_mpq_lar_core_solver.m_d_basis = m_mpq_lar_core_solver.m_r_basis;
+
+	/*
+    unsigned j = A_d().column_count();
+    A_d().add_column();
+    lean_assert(m_mpq_lar_core_solver.m_d_x.size() == j);
+    //        lean_assert(m_mpq_lar_core_solver.m_d_low_bounds.size() == j && m_mpq_lar_core_solver.m_d_upper_bounds.size() == j);  // restore later
+    m_mpq_lar_core_solver.m_d_x.resize(j + 1 ); 
+    m_mpq_lar_core_solver.m_d_low_bounds.resize(j + 1);
+    m_mpq_lar_core_solver.m_d_upper_bounds.resize(j + 1);
+    lean_assert(m_mpq_lar_core_solver.m_d_heading.size() == j); // as A().column_count() on the entry to the method
+    if (register_in_basis) {
+        A_d().add_row();
+        m_mpq_lar_core_solver.m_d_heading.push_back(m_mpq_lar_core_solver.m_d_basis.size());
+        m_mpq_lar_core_solver.m_d_basis.push_back(j);
+    }else {
+        m_mpq_lar_core_solver.m_d_heading.push_back(- static_cast<int>(m_mpq_lar_core_solver.m_d_nbasis.size()) - 1);
+        m_mpq_lar_core_solver.m_d_nbasis.push_back(j);
+        }*/
+}
+
+void adjust_initial_state_for_tableau_rows() {
+    for (unsigned j = 0; j < m_terms.size(); j++) {
+        if (contains(m_ext_vars_to_columns, j + m_terms_start_index))
+            continue;
+        add_row_from_term_no_constraint(m_terms[j], j + m_terms_start_index);
+    }
+}
+
+// this fills the last row of A_d and sets the basis column: -1 in the last column of the row
+void fill_last_row_of_A_d(static_matrix<double, double> & A, const lar_term* ls) {
+    lean_assert(A.row_count() > 0);
+    lean_assert(A.column_count() > 0);
+    unsigned last_row = A.row_count() - 1;
+    lean_assert(A.m_rows[last_row].empty());
+
+    for (auto & t : ls->m_coeffs) {
+        lean_assert(!is_zero(t.second));
+        var_index j = t.first;
+        A.set(last_row, j, - t.second.get_double());
+    }
+
+    unsigned basis_j = A.column_count() - 1;
+    A.set(last_row, basis_j, - 1 );
+}
+
+void update_free_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index constr_ind) {
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::upper_bound;
+        lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
+        lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+        }
+        set_upper_bound_witness(j, constr_ind);
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::low_bound;
+        lean_assert(m_mpq_lar_core_solver.m_r_upper_bounds.size() > j);
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+        }
+        set_low_bound_witness(j, constr_ind);
+        break;
+    case EQ:
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+        m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+        set_upper_bound_witness(j, constr_ind);
+        set_low_bound_witness(j, constr_ind);
+        break;
+
+    default:
+        lean_unreachable();
+                
+    }
+    m_columns_with_changed_bound.insert(j);
+}
+
+void update_upper_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::upper_bound);
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_upper_bounds()[j]) {
+                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+                set_upper_bound_witness(j, ci);
+                m_columns_with_changed_bound.insert(j);
+            }
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        m_mpq_lar_core_solver.m_column_types[j] = column_type::boxed;
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+            set_low_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
+            }                     
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                set_low_bound_witness(j, ci);
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+            break;
+        }
+        break;
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+    
+void update_boxed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::boxed && m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            if (up < m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+                set_upper_bound_witness(j, ci);
+                m_columns_with_changed_bound.insert(j);
+            }
+
+            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                lean_assert(false);
+                m_infeasible_column_index = j;
+            } else {
+                if (m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j])
+                    m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }                    
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+            }
+            if (low > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else if ( low == m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);                    
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+                m_columns_with_changed_bound.insert(j);
+            }
+                
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+void update_low_bound_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_mpq_lar_core_solver.m_column_types()[j] == column_type::low_bound);
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        y_of_bound = -1;
+    case LE:
+        {
+            auto up = numeric_pair<mpq>(right_side, y_of_bound);
+            m_mpq_lar_core_solver.m_r_upper_bounds[j] = up;
+            set_upper_bound_witness(j, ci);
+            m_columns_with_changed_bound.insert(j);
+
+            if (up < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+            } else {
+                m_mpq_lar_core_solver.m_column_types[j] = m_mpq_lar_core_solver.m_r_low_bounds()[j] < m_mpq_lar_core_solver.m_r_upper_bounds()[j]? column_type::boxed : column_type::fixed;
+            }                    
+        }
+        break;
+    case GT:
+        y_of_bound = 1;
+    case GE:            
+        {
+            auto low = numeric_pair<mpq>(right_side, y_of_bound);
+            if (low > m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = low;
+                m_columns_with_changed_bound.insert(j);
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case EQ:
+        {
+            auto v = numeric_pair<mpq>(right_side, zero_of_type<mpq>());
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else {
+                m_mpq_lar_core_solver.m_r_low_bounds[j] = m_mpq_lar_core_solver.m_r_upper_bounds[j] = v;
+                set_low_bound_witness(j, ci);
+                set_upper_bound_witness(j, ci);
+                m_mpq_lar_core_solver.m_column_types[j] = column_type::fixed;
+            }
+            m_columns_with_changed_bound.insert(j);
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+
+void update_fixed_column_type_and_bound(var_index j, lconstraint_kind kind, const mpq & right_side, constraint_index ci) {
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_column_types()[j] == column_type::fixed && m_mpq_lar_core_solver.m_r_low_bounds()[j] == m_mpq_lar_core_solver.m_r_upper_bounds()[j]));
+    lean_assert(m_status == INFEASIBLE || (m_mpq_lar_core_solver.m_r_low_bounds()[j].y.is_zero() && m_mpq_lar_core_solver.m_r_upper_bounds()[j].y.is_zero()));
+    auto v = numeric_pair<mpq>(right_side, mpq(0));
+        
+    mpq y_of_bound(0);
+    switch (kind) {
+    case LT:
+        if (v <= m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+            m_status = INFEASIBLE;
+            m_infeasible_column_index = j;
+            set_upper_bound_witness(j, ci);
+        }                   
+        break;
+    case LE:
+        {
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);
+            }                   
+        }
+        break;
+    case GT:
+        {
+            if (v >= m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index =j;
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case GE:            
+        {
+            if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);
+            }
+        }
+        break;
+    case EQ:
+        {
+            if (v < m_mpq_lar_core_solver.m_r_low_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_upper_bound_witness(j, ci);                    
+            } else if (v > m_mpq_lar_core_solver.m_r_upper_bounds[j]) {
+                m_status = INFEASIBLE;
+                m_infeasible_column_index = j;
+                set_low_bound_witness(j, ci);                    
+            } 
+            break;
+        }
+
+    default:
+        lean_unreachable();
+                
+    }
+}
+    
diff --git a/src/util/lp/iterator_on_column.h b/src/util/lp/iterator_on_column.h
index 5ad925d3e..215514b39 100644
--- a/src/util/lp/iterator_on_column.h
+++ b/src/util/lp/iterator_on_column.h
@@ -11,7 +11,7 @@ template <typename T, typename X>
 struct iterator_on_column:linear_combination_iterator<T> {
     const vector<column_cell>& m_column; // the offset in term coeffs
     const static_matrix<T, X> & m_A;
-    int m_i = -1; // the initial offset in the column
+    int m_i; // the initial offset in the column
     unsigned size() const { return m_column.size(); }
     iterator_on_column(const vector<column_cell>& column, const static_matrix<T,X> & A) // the offset in term coeffs
         :
diff --git a/src/util/lp/iterator_on_indexed_vector.h b/src/util/lp/iterator_on_indexed_vector.h
index c19c08e64..532b62617 100644
--- a/src/util/lp/iterator_on_indexed_vector.h
+++ b/src/util/lp/iterator_on_indexed_vector.h
@@ -8,8 +8,11 @@ namespace lean {
 template <typename T>
 struct iterator_on_indexed_vector:linear_combination_iterator<T> {
     const indexed_vector<T> & m_v;
-    unsigned m_offset = 0;
-    iterator_on_indexed_vector(const indexed_vector<T> & v) : m_v(v){}
+    unsigned m_offset;
+    iterator_on_indexed_vector(const indexed_vector<T> & v) :
+        m_v(v),
+        m_offset(0)
+    {}
     unsigned size() const { return m_v.m_index.size(); }
     bool next(T & a, unsigned & i) {
         if (m_offset >= m_v.m_index.size())
diff --git a/src/util/lp/iterator_on_pivot_row.h b/src/util/lp/iterator_on_pivot_row.h
index 19c0c7df1..1a9381a70 100644
--- a/src/util/lp/iterator_on_pivot_row.h
+++ b/src/util/lp/iterator_on_pivot_row.h
@@ -7,12 +7,14 @@
 namespace lean {
 template <typename T>
 struct iterator_on_pivot_row:linear_combination_iterator<T> {
-    bool m_basis_returned = false;
+    bool m_basis_returned;
     const indexed_vector<T> & m_v;
     unsigned m_basis_j;
     iterator_on_indexed_vector<T> m_it;
     unsigned size() const { return m_it.size(); }
-    iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) : m_v(v), m_basis_j(basis_j), m_it(v) {}
+    iterator_on_pivot_row(const indexed_vector<T> & v, unsigned basis_j) :
+        m_basis_returned(false),
+        m_v(v), m_basis_j(basis_j), m_it(v) {}
     bool next(T & a, unsigned & i) {
         if (m_basis_returned == false) {
             m_basis_returned = true;
diff --git a/src/util/lp/iterator_on_row.h b/src/util/lp/iterator_on_row.h
index 212768db9..96a1a8cf3 100644
--- a/src/util/lp/iterator_on_row.h
+++ b/src/util/lp/iterator_on_row.h
@@ -8,8 +8,8 @@ namespace lean {
 template <typename T>
 struct iterator_on_row:linear_combination_iterator<T> {
     const vector<row_cell<T>> & m_row;
-    unsigned m_i= 0; // offset
-    iterator_on_row(const vector<row_cell<T>> & row) : m_row(row)
+    unsigned m_i; // offset
+    iterator_on_row(const vector<row_cell<T>> & row) : m_row(row), m_i(0)
     {}
     unsigned size() const { return m_row.size(); }
     bool next(T & a, unsigned & i) {
diff --git a/src/util/lp/iterator_on_term_with_basis_var.h b/src/util/lp/iterator_on_term_with_basis_var.h
index 9279dd637..74a2c4038 100644
--- a/src/util/lp/iterator_on_term_with_basis_var.h
+++ b/src/util/lp/iterator_on_term_with_basis_var.h
@@ -9,11 +9,12 @@
 namespace lean {
 struct iterator_on_term_with_basis_var:linear_combination_iterator<mpq> {
     std::unordered_map<unsigned, mpq>::const_iterator m_i; // the offset in term coeffs
-    bool m_term_j_returned = false;
+    bool m_term_j_returned;
     const lar_term & m_term;
     unsigned m_term_j;
     unsigned size() const {return static_cast<unsigned>(m_term.m_coeffs.size() + 1);}
     iterator_on_term_with_basis_var(const lar_term & t, unsigned term_j) :
+        m_term_j_returned(false),
         m_i(t.m_coeffs.begin()),
         m_term(t),
         m_term_j(term_j) {}
diff --git a/src/util/lp/lar_core_solver.h b/src/util/lp/lar_core_solver.h
index 00abb188d..58151c8be 100644
--- a/src/util/lp/lar_core_solver.h
+++ b/src/util/lp/lar_core_solver.h
@@ -25,7 +25,7 @@ class lar_core_solver  {
     // to grow and is set to -1  otherwise
     int m_sign_of_entering_delta;
     vector<std::pair<mpq, unsigned>> m_infeasible_linear_combination;
-    int m_infeasible_sum_sign = 0; // todo: get rid of this field
+    int m_infeasible_sum_sign; // todo: get rid of this field
     vector<numeric_pair<mpq>> m_right_sides_dummy;
     vector<mpq> m_costs_dummy;
     vector<double> m_d_right_sides_dummy;
diff --git a/src/util/lp/lar_core_solver.hpp b/src/util/lp/lar_core_solver.hpp
index ded6762c5..a6dd7e3e0 100644
--- a/src/util/lp/lar_core_solver.hpp
+++ b/src/util/lp/lar_core_solver.hpp
@@ -14,7 +14,8 @@ namespace lean {
 lar_core_solver::lar_core_solver(
                                  lp_settings & settings,
                                  const column_namer & column_names
-):
+                                 ):
+    m_infeasible_sum_sign(0),
     m_r_solver(m_r_A,
                     m_right_sides_dummy,
                     m_r_x,
diff --git a/src/util/lp/lar_solver.h b/src/util/lp/lar_solver.h
index 39961c3ad..b8609675a 100644
--- a/src/util/lp/lar_solver.h
+++ b/src/util/lp/lar_solver.h
@@ -35,7 +35,7 @@ namespace lean {
 class lar_solver : public column_namer {
     //////////////////// fields //////////////////////////
     lp_settings m_settings;
-    stacked_value<lp_status> m_status = OPTIMAL;
+    stacked_value<lp_status> m_status;
     stacked_value<simplex_strategy_enum> m_simplex_strategy;
 	std::unordered_map<unsigned, var_index> m_ext_vars_to_columns;
     vector<unsigned> m_columns_to_ext_vars_or_term_indices;
@@ -46,13 +46,13 @@ class lar_solver : public column_namer {
     int_set m_columns_with_changed_bound;
     int_set m_rows_with_changed_bounds;
 	int_set m_basic_columns_with_changed_cost;
-    stacked_value<int> m_infeasible_column_index = -1; // such can be found at the initialization step
+    stacked_value<int> m_infeasible_column_index; // such can be found at the initialization step
     stacked_value<unsigned> m_term_count;
     vector<lar_term*> m_terms;
     vector<lar_term*> m_orig_terms;
-    const var_index m_terms_start_index = 1000000;
+    const var_index m_terms_start_index;
     indexed_vector<mpq> m_column_buffer;    
-    std::function<column_type (unsigned)> m_column_type_function = [this] (unsigned j) {return m_mpq_lar_core_solver.m_column_types()[j];};    
+    std::function<column_type (unsigned)> m_column_type_function;    
 public:
     lar_core_solver m_mpq_lar_core_solver;
 	unsigned constraint_count() const {
@@ -83,7 +83,12 @@ public:
     lar_solver() : m_mpq_lar_core_solver(
                                          m_settings,
                                          *this
-                                         ) {}
+                                         ),
+                   m_status(OPTIMAL),
+                   m_infeasible_column_index(-1),
+                   m_terms_start_index(1000000),
+                   m_column_type_function ([this] (unsigned j) {return m_mpq_lar_core_solver.m_column_types()[j];})    
+    {}
     
     void set_propagate_bounds_on_pivoted_rows_mode(bool v) {
         m_mpq_lar_core_solver.m_r_solver.m_pivoted_rows = v? (& m_rows_with_changed_bounds) : nullptr;
diff --git a/src/util/lp/linear_combination_iterator.h b/src/util/lp/linear_combination_iterator.h
index 9a1a2de77..634accfd4 100644
--- a/src/util/lp/linear_combination_iterator.h
+++ b/src/util/lp/linear_combination_iterator.h
@@ -16,7 +16,7 @@ struct linear_combination_iterator {
 template <typename T>
 struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
     vector<std::pair<T, unsigned>> & m_vector;
-    int m_offset = 0;
+    int m_offset;
     bool next(T & a, unsigned & i) {
         if(m_offset >= m_vector.size())
             return false;
@@ -40,7 +40,10 @@ struct linear_combination_iterator_on_vector : linear_combination_iterator<T> {
     linear_combination_iterator<T> * clone() {
         return new linear_combination_iterator_on_vector(m_vector);
     }
-    linear_combination_iterator_on_vector(vector<std::pair<T, unsigned>> & vec): m_vector(vec) {}
+    linear_combination_iterator_on_vector(vector<std::pair<T, unsigned>> & vec):
+        m_vector(vec),
+        m_offset(0)
+    {}
     unsigned size() const { return m_vector.size(); }
 };
 
diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h
index d8618f32a..86382dd49 100644
--- a/src/util/lp/lp_core_solver_base.h
+++ b/src/util/lp/lp_core_solver_base.h
@@ -17,7 +17,7 @@ namespace lean {
 
 template <typename T, typename X> // X represents the type of the x variable and the bounds
 class lp_core_solver_base {    
-    unsigned m_total_iterations = 0;
+    unsigned m_total_iterations;
     unsigned inc_total_iterations() { ++m_settings.st().m_total_iterations; return m_total_iterations++; }
 private:
     lp_status m_status;
@@ -25,7 +25,7 @@ public:
     bool current_x_is_feasible() const { return m_inf_set.size() == 0; }
     bool current_x_is_infeasible() const { return m_inf_set.size() != 0; }
     int_set m_inf_set;
-    bool m_using_infeas_costs = false;
+    bool m_using_infeas_costs;
 
 
     vector<unsigned> m_columns_nz; // m_columns_nz[i] keeps an approximate value of non zeroes the i-th column
@@ -42,23 +42,23 @@ public:
     lp_settings & m_settings;
     vector<T> m_y; // the buffer for yB = cb
     // a device that is able to solve Bx=c, xB=d, and change the basis
-    lu<T, X> * m_factorization = nullptr;
+    lu<T, X> * m_factorization;
     const column_namer & m_column_names;
     indexed_vector<T> m_w; // the vector featuring in 24.3 of the Chvatal book
     vector<T> m_d; // the vector of reduced costs
     indexed_vector<T> m_ed; // the solution of B*m_ed = a
-    unsigned m_iters_with_no_cost_growing = 0;
+    unsigned m_iters_with_no_cost_growing;
     const vector<column_type> & m_column_types;
     const vector<X> & m_low_bounds;
     const vector<X> & m_upper_bounds;
     vector<T> m_column_norms; // the approximate squares of column norms that help choosing a profitable column
     vector<X> m_copy_of_xB;
-    unsigned m_basis_sort_counter = 0;
+    unsigned m_basis_sort_counter;
     vector<T> m_steepest_edge_coefficients;
     vector<unsigned> m_trace_of_basis_change_vector; // the even positions are entering, the odd positions are leaving
-    bool m_tracing_basis_changes = false;
-    int_set* m_pivoted_rows = nullptr;
-    bool m_look_for_feasible_solution_only = false;
+    bool m_tracing_basis_changes;
+    int_set* m_pivoted_rows;
+    bool m_look_for_feasible_solution_only;
     void start_tracing_basis_changes() {
         m_trace_of_basis_change_vector.resize(0);
         m_tracing_basis_changes = true;
diff --git a/src/util/lp/lp_core_solver_base.hpp b/src/util/lp/lp_core_solver_base.hpp
index 4b67e86ee..5ae1e399c 100644
--- a/src/util/lp/lp_core_solver_base.hpp
+++ b/src/util/lp/lp_core_solver_base.hpp
@@ -45,7 +45,14 @@ lp_core_solver_base(static_matrix<T, X> & A,
     m_upper_bounds(upper_bound_values),
     m_column_norms(m_n()),
     m_copy_of_xB(m_m()),
-    m_steepest_edge_coefficients(A.column_count()) {
+    m_steepest_edge_coefficients(A.column_count()),
+    m_total_iterations(0),
+    m_using_infeas_costs(false),
+    m_iters_with_no_cost_growing(0),
+    m_basis_sort_counter(0),
+    m_tracing_basis_changes(false),
+    m_pivoted_rows(nullptr),
+    m_look_for_feasible_solution_only(false) {
     lean_assert(bounds_for_boxed_are_set_correctly());    
     init();
     init_basis_heading_and_non_basic_columns_vector();
diff --git a/src/util/lp/lp_dual_simplex.h b/src/util/lp/lp_dual_simplex.h
index 78ff08b0a..4dff2a4f1 100644
--- a/src/util/lp/lp_dual_simplex.h
+++ b/src/util/lp/lp_dual_simplex.h
@@ -11,7 +11,7 @@ namespace lean {
 
 template <typename T, typename X>
 class lp_dual_simplex: public lp_solver<T, X> {
-    lp_dual_core_solver<T, X> * m_core_solver = nullptr;
+    lp_dual_core_solver<T, X> * m_core_solver;
     vector<T> m_b_copy;
     vector<T> m_low_bounds; // We don't have a convention here that all low bounds are zeros. At least it does not hold for the first stage solver
     vector<column_type> m_column_types_of_core_solver;
@@ -24,6 +24,7 @@ public:
         }
     }
 
+    lp_dual_simplex() : m_core_solver(nullptr) {}
 
 
     void decide_on_status_after_stage1();
diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h
index 62464ca8b..e4d869289 100644
--- a/src/util/lp/lp_primal_core_solver.h
+++ b/src/util/lp/lp_primal_core_solver.h
@@ -38,13 +38,13 @@ public:
     vector<breakpoint<X>> m_breakpoints;
     binary_heap_priority_queue<X> m_breakpoint_indices_queue;
     indexed_vector<T> m_beta; // see Swietanowski working vector beta for column norms
-    T m_epsilon_of_reduced_cost = T(1)/T(10000000);
+    T m_epsilon_of_reduced_cost;
     vector<T> m_costs_backup;
     T m_converted_harris_eps;
     unsigned m_inf_row_index_for_tableau;
     bool m_bland_mode_tableau;
     int_set m_left_basis_tableau;
-    unsigned m_bland_mode_threshold = 1000;
+    unsigned m_bland_mode_threshold;
     unsigned m_left_basis_repeated;
     vector<unsigned> m_leaving_candidates;
     //    T m_converted_harris_eps = convert_struct<T, double>::convert(this->m_settings.harris_feasibility_tolerance);
@@ -905,6 +905,8 @@ public:
                                   column_type_array,
                                   low_bound_values,
                                   upper_bound_values),
+        m_epsilon_of_reduced_cost(T(1)/T(10000000)),
+        m_bland_mode_threshold(1000),
         m_beta(A.row_count()) {
 
         if (!(numeric_traits<T>::precise())) {
diff --git a/src/util/lp/lp_primal_simplex.h b/src/util/lp/lp_primal_simplex.h
index 3b1288fc7..fcddb4eb1 100644
--- a/src/util/lp/lp_primal_simplex.h
+++ b/src/util/lp/lp_primal_simplex.h
@@ -15,7 +15,7 @@
 namespace lean {
 template <typename T, typename X>
 class lp_primal_simplex: public lp_solver<T, X> {
-    lp_primal_core_solver<T, X> * m_core_solver = nullptr;
+    lp_primal_core_solver<T, X> * m_core_solver;
     vector<X> m_low_bounds;
 private:
     unsigned original_rows() { return this->m_external_rows_to_core_solver_rows.size(); }
@@ -28,7 +28,7 @@ private:
 
     void set_scaled_costs();
 public:
-    lp_primal_simplex() {}
+    lp_primal_simplex(): m_core_solver(nullptr) {}
 
     column_info<T> * get_or_create_column_info(unsigned column);
 
diff --git a/src/util/lp/lp_settings.h b/src/util/lp/lp_settings.h
index 5fab4e451..e97f5f5f5 100644
--- a/src/util/lp/lp_settings.h
+++ b/src/util/lp/lp_settings.h
@@ -106,49 +106,49 @@ private:
     default_lp_resource_limit m_default_resource_limit;
     lp_resource_limit* m_resource_limit;
     // used for debug output
-    std::ostream* m_debug_out = &std::cout;
+    std::ostream* m_debug_out;
     // used for messages, for example, the computation progress messages
-    std::ostream* m_message_out = &std::cout;
+    std::ostream* m_message_out;
 
     stats  m_stats;
 
 public:
-    unsigned reps_in_scaler = 20;
+    unsigned reps_in_scaler;
     // when the absolute value of an element is less than pivot_epsilon
     // in pivoting, we treat it as a zero
-    double pivot_epsilon = 0.00000001;
+    double pivot_epsilon;
     // see Chatal, page 115
-    double positive_price_epsilon = 1e-7;
+    double positive_price_epsilon;
     // a quatation "if some choice of the entering vairable leads to an eta matrix
     // whose diagonal element in the eta column is less than e2 (entering_diag_epsilon) in magnitude, the this choice is rejected ...
-    double entering_diag_epsilon = 1e-8;
-    int c_partial_pivoting = 10; // this is the constant c from page 410
-    unsigned depth_of_rook_search = 4;
-    bool using_partial_pivoting = true;
+    double entering_diag_epsilon;
+    int c_partial_pivoting; // this is the constant c from page 410
+    unsigned depth_of_rook_search;
+    bool using_partial_pivoting;
     // dissertation of Achim Koberstein
     // if Bx - b is different at any component more that refactor_epsilon then we refactor
-    double refactor_tolerance = 1e-4;
-    double pivot_tolerance = 1e-6;
-    double zero_tolerance = 1e-12;
-    double drop_tolerance = 1e-14;
-    double tolerance_for_artificials = 1e-4;
-    double can_be_taken_to_basis_tolerance = 0.00001;
+    double refactor_tolerance;
+    double pivot_tolerance;
+    double zero_tolerance;
+    double drop_tolerance;
+    double tolerance_for_artificials;
+    double can_be_taken_to_basis_tolerance;
 
-    unsigned percent_of_entering_to_check = 5; // we try to find a profitable column in a percentage of the columns
-    bool use_scaling = true;
-    double scaling_maximum = 1;
-    double scaling_minimum = 0.5;
-    double harris_feasibility_tolerance = 1e-7; // page 179 of Istvan Maros
-    double ignore_epsilon_of_harris = 10e-5;
-    unsigned max_number_of_iterations_with_no_improvements = 2000000;
-    unsigned max_total_number_of_iterations = 20000000;
-    double time_limit = std::numeric_limits<double>::max(); // the maximum time limit of the total run time in seconds
+    unsigned percent_of_entering_to_check; // we try to find a profitable column in a percentage of the columns
+    bool use_scaling;
+    double scaling_maximum;
+    double scaling_minimum;
+    double harris_feasibility_tolerance; // page 179 of Istvan Maros
+    double ignore_epsilon_of_harris;
+    unsigned max_number_of_iterations_with_no_improvements;
+    unsigned max_total_number_of_iterations;
+    double time_limit; // the maximum time limit of the total run time in seconds
     // dual section
-    double dual_feasibility_tolerance = 1e-7; // // page 71 of the PhD thesis of Achim Koberstein
-    double primal_feasibility_tolerance = 1e-7; // page 71 of the PhD thesis of Achim Koberstein
-    double relative_primal_feasibility_tolerance = 1e-9; // page 71 of the PhD thesis of Achim Koberstein
+    double dual_feasibility_tolerance; // // page 71 of the PhD thesis of Achim Koberstein
+    double primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
+    double relative_primal_feasibility_tolerance; // page 71 of the PhD thesis of Achim Koberstein
 
-    bool m_bound_propagation = true;
+    bool m_bound_propagation;
     
     bool bound_progation() const {
         return m_bound_propagation;
@@ -158,7 +158,53 @@ public:
         return m_bound_propagation;
     }
     
-    lp_settings() : m_default_resource_limit(*this), m_resource_limit(&m_default_resource_limit) {}
+    lp_settings() : m_default_resource_limit(*this),
+                    m_resource_limit(&m_default_resource_limit),
+                    m_debug_out( &std::cout),
+                    m_message_out(&std::cout),
+                    reps_in_scaler(20),
+                    pivot_epsilon(0.00000001),
+                    positive_price_epsilon(1e-7),
+                    entering_diag_epsilon ( 1e-8),
+                    c_partial_pivoting ( 10), // this is the constant c from page 410
+                    depth_of_rook_search ( 4),
+                    using_partial_pivoting ( true),
+                    // dissertation of Achim Koberstein
+                    // if Bx - b is different at any component more that refactor_epsilon then we refactor
+                    refactor_tolerance ( 1e-4),
+                    pivot_tolerance ( 1e-6),
+                    zero_tolerance ( 1e-12),
+                    drop_tolerance ( 1e-14),
+                    tolerance_for_artificials ( 1e-4),
+                    can_be_taken_to_basis_tolerance ( 0.00001),
+                  
+                    percent_of_entering_to_check ( 5),// we try to find a profitable column in a percentage of the columns
+                    use_scaling ( true),
+                    scaling_maximum ( 1),
+                    scaling_minimum ( 0.5),
+                    harris_feasibility_tolerance ( 1e-7), // page 179 of Istvan Maros
+                    ignore_epsilon_of_harris ( 10e-5),
+                    max_number_of_iterations_with_no_improvements ( 2000000),
+                    max_total_number_of_iterations ( 20000000),
+                    time_limit ( std::numeric_limits<double>::max()), // the maximum time limit of the total run time in seconds
+                    // dual section
+                    dual_feasibility_tolerance ( 1e-7), // // page 71 of the PhD thesis of Achim Koberstein
+                    primal_feasibility_tolerance ( 1e-7), // page 71 of the PhD thesis of Achim Koberstein
+                    relative_primal_feasibility_tolerance ( 1e-9), // page 71 of the PhD thesis of Achim Koberstein
+                    m_bound_propagation ( true),
+                    presolve_with_double_solver_for_lar(true),
+                    m_simplex_strategy(simplex_strategy_enum::tableau_rows),
+                  report_frequency(1000),
+                  print_statistics(false),
+                  column_norms_update_frequency(12000),
+                  scale_with_ratio(true),
+                  density_threshold(0.7),
+                  use_breakpoints_in_feasibility_search(false),
+                  random_seed(1),
+                  max_row_length_for_bound_propagation(300),
+                  backup_costs(true),
+                  column_number_threshold_for_using_lu_in_lar_solver(4000)
+    {}
 
     void set_resource_limit(lp_resource_limit& lim) { m_resource_limit = &lim; }
     bool get_cancel_flag() const { return m_resource_limit->get_cancel_flag(); }
@@ -227,8 +273,8 @@ public:
         return is_eps_small_general<T>(t, tolerance_for_artificials);
     }
     // the method of lar solver to use
-    bool presolve_with_double_solver_for_lar = true;
-    simplex_strategy_enum m_simplex_strategy = simplex_strategy_enum::tableau_rows;
+    bool presolve_with_double_solver_for_lar;
+    simplex_strategy_enum m_simplex_strategy;
     simplex_strategy_enum simplex_strategy() const {
         return m_simplex_strategy;
     }
@@ -250,20 +296,20 @@ public:
         return m_simplex_strategy == simplex_strategy_enum::tableau_rows;
     }
     
-    int report_frequency = 1000;
-    bool print_statistics = false;
-    unsigned column_norms_update_frequency = 12000;
-    bool scale_with_ratio = true;
-    double density_threshold = 0.7; // need to tune it up, todo
+    int report_frequency;
+    bool print_statistics;
+    unsigned column_norms_update_frequency;
+    bool scale_with_ratio;
+    double density_threshold; // need to tune it up, todo
 #ifdef LEAN_DEBUG
     static unsigned ddd; // used for debugging    
 #endif
-    bool use_breakpoints_in_feasibility_search = false;
-    unsigned random_seed = 1;
+    bool use_breakpoints_in_feasibility_search;
+    unsigned random_seed;
     static unsigned long random_next;
-    unsigned max_row_length_for_bound_propagation = 300;
-    bool backup_costs = true;
-    unsigned column_number_threshold_for_using_lu_in_lar_solver = 4000;
+    unsigned max_row_length_for_bound_propagation;
+    bool backup_costs;
+    unsigned column_number_threshold_for_using_lu_in_lar_solver;
 }; // end of lp_settings class
 
 
diff --git a/src/util/lp/lp_solver.h b/src/util/lp/lp_solver.h
index eeb3ff6d3..1bfe7dcdc 100644
--- a/src/util/lp/lp_solver.h
+++ b/src/util/lp/lp_solver.h
@@ -39,10 +39,10 @@ protected:
     T get_column_cost_value(unsigned j, column_info<T> * ci) const;
 public:
     unsigned m_total_iterations;
-    static_matrix<T, X>* m_A = nullptr; // this is the matrix of constraints
+    static_matrix<T, X>* m_A; // this is the matrix of constraints
     vector<T> m_b; // the right side vector
-    unsigned m_first_stage_iterations = 0;
-    unsigned m_second_stage_iterations = 0;
+    unsigned m_first_stage_iterations;
+    unsigned m_second_stage_iterations;
     std::unordered_map<unsigned, lp_constraint<T, X>> m_constraints;
     std::unordered_map<var_index, column_info<T>*> m_map_from_var_index_to_column_info;
     std::unordered_map<unsigned, std::unordered_map<unsigned, T> > m_A_values;
@@ -52,8 +52,8 @@ public:
     std::unordered_map<unsigned, unsigned> m_core_solver_columns_to_external_columns;
     vector<T> m_column_scale;
     std::unordered_map<unsigned, std::string>  m_name_map;
-    unsigned m_artificials = 0;
-    unsigned m_slacks = 0;
+    unsigned m_artificials;
+    unsigned m_slacks;
     vector<column_type> m_column_types;
     vector<T> m_costs;
     vector<T> m_x;
@@ -63,10 +63,17 @@ public:
     vector<int> m_heading;
 
 
-    lp_status m_status = lp_status::UNKNOWN;
+    lp_status m_status;
 
     lp_settings m_settings;
-    lp_solver() {}
+    lp_solver():
+        m_A(nullptr), // this is the matrix of constraints
+        m_first_stage_iterations (0),
+        m_second_stage_iterations (0),
+        m_artificials (0),
+        m_slacks (0),
+        m_status(lp_status::UNKNOWN)
+    {}
 
     unsigned row_count() const { return this->m_A->row_count(); }
 
@@ -232,14 +239,6 @@ protected:
         out << "extended A[" << this->m_A->row_count() << "," << this->m_A->column_count() << "]" << std::endl;
     }
 
-    struct row_tighten_stats {
-        unsigned n_of_new_bounds = 0;
-        unsigned n_of_fixed = 0;
-        bool is_obsolete = false;
-    };
-    
-
-
 public:
     lp_settings & settings() { return m_settings;}
     void print_model(std::ostream & s) const {
diff --git a/src/util/lp/lu.h b/src/util/lp/lu.h
index 415b7f978..f58c0d845 100644
--- a/src/util/lp/lu.h
+++ b/src/util/lp/lu.h
@@ -110,7 +110,7 @@ enum class LU_status { OK, Degenerated};
 // Using Suhl-Suhl method described in the dissertation of Achim Koberstein, Chapter 5
 template <typename T, typename X>
 class lu {
-    LU_status m_status = LU_status::OK;
+    LU_status m_status;
 public:
     // the fields
     unsigned m_dim;
@@ -123,12 +123,12 @@ public:
     
     vector<tail_matrix<T, X> *> m_tail;
     lp_settings & m_settings;
-    bool m_failure = false;
+    bool m_failure;
     indexed_vector<T> m_row_eta_work_vector;
     indexed_vector<T> m_w_for_extension;
     indexed_vector<T> m_y_copy;
     indexed_vector<unsigned> m_ii; //to optimize the work with the m_index fields
-    unsigned m_refactor_counter = 0;
+    unsigned m_refactor_counter;
     // constructor
     // if A is an m by n matrix then basis has length m and values in [0,n); the values are all different
     // they represent the set of m columns
diff --git a/src/util/lp/lu.hpp b/src/util/lp/lu.hpp
index ba5e092d3..c0b3d19ad 100644
--- a/src/util/lp/lu.hpp
+++ b/src/util/lp/lu.hpp
@@ -118,7 +118,10 @@ lu<T, X>::lu(static_matrix<T, X> const & A,
     m_r_wave(m_dim),
     m_U(A, basis), // create the square matrix that eventually will be factorized
     m_settings(settings),
-    m_row_eta_work_vector(A.row_count()){
+    m_row_eta_work_vector(A.row_count()),
+    m_status(LU_status::OK),
+    m_failure(false),
+    m_refactor_counter(0) {
     lean_assert(!(numeric_traits<T>::precise() && settings.use_tableau()));
 #ifdef LEAN_DEBUG
     debug_test_of_basis(A, basis);
diff --git a/src/util/lp/mps_reader.h b/src/util/lp/mps_reader.h
index c79a7d426..206220584 100644
--- a/src/util/lp/mps_reader.h
+++ b/src/util/lp/mps_reader.h
@@ -93,22 +93,28 @@ template <typename T, typename X>
 class mps_reader {
     enum row_type { Cost, Less_or_equal, Greater_or_equal, Equal };
     struct bound {
-        bool m_low_is_set = true;
+        bool m_low_is_set;
         T m_low;
-        bool m_upper_is_set = false;
+        bool m_upper_is_set;
         T m_upper;
-        bool m_value_is_fixed = false;
+        bool m_value_is_fixed;
         T m_fixed_value;
-        bool m_free = false;
+        bool m_free;
         // constructor
-        bound() : m_low(numeric_traits<T>::zero()) {} // it seems all mps files I have seen have the default low value 0 on a variable
+        bound() : m_low(numeric_traits<T>::zero()),
+                  m_low_is_set(true),
+                  m_upper_is_set(false),
+                  m_value_is_fixed(false),
+                  m_free(false) {} // it seems all mps files I have seen have the default low value 0 on a variable
     };
 
     struct column {
         std::string m_name;
-        bound * m_bound = nullptr;
+        bound * m_bound;
         unsigned m_index;
-        column(std::string name, unsigned index): m_name(name), m_index(index) {
+        column(std::string name, unsigned index): m_name(name),
+                                                  m_bound(nullptr),
+                                                  m_index(index) {
         }
     };
 
@@ -116,15 +122,18 @@ class mps_reader {
         row_type m_type;
         std::string m_name;
         std::unordered_map<std::string, T> m_row_columns;
-        T m_right_side = numeric_traits<T>::zero();
+        T m_right_side;
         unsigned m_index;
-        T m_range = numeric_traits<T>::zero();
-        row(row_type type, std::string name, unsigned index) : m_type(type), m_name(name), m_index(index) {
+        T m_range;
+        row(row_type type, std::string name, unsigned index) : m_type(type), m_name(name), m_index(index),
+                                                               m_right_side(zero_of_type<T>()),
+                                                               m_range(zero_of_type<T>())
+        {
         }
     };
 
     std::string m_file_name;
-    bool m_is_OK = true;
+    bool m_is_OK;
     std::unordered_map<std::string, row *> m_rows;
     std::unordered_map<std::string, column *> m_columns;
     std::unordered_map<std::string, unsigned> m_names_to_var_index;
@@ -133,9 +142,9 @@ class mps_reader {
     std::string m_cost_row_name;
     std::ifstream m_file_stream;
     // needed to adjust the index row
-    unsigned m_cost_line_count = 0;
-    unsigned m_line_number = 0;
-    std::ostream * m_message_stream = & std::cout;
+    unsigned m_cost_line_count;
+    unsigned m_line_number;
+    std::ostream * m_message_stream;
 
     void set_m_ok_to_false() {
         *m_message_stream << "setting m_is_OK to false" << std::endl;
@@ -737,8 +746,11 @@ public:
     }
 
     mps_reader(std::string file_name):
-        m_file_name(file_name), m_file_stream(file_name) {
-    }
+        m_is_OK(true),
+        m_file_name(file_name), m_file_stream(file_name),
+        m_cost_line_count(0),
+        m_line_number(0),
+        m_message_stream(& std::cout) {}
     void read() {
         if (!m_file_stream.is_open()){
             set_m_ok_to_false();
diff --git a/src/util/lp/numeric_pair.h b/src/util/lp/numeric_pair.h
index 2cea4158d..84c99b3b1 100644
--- a/src/util/lp/numeric_pair.h
+++ b/src/util/lp/numeric_pair.h
@@ -101,16 +101,14 @@ struct numeric_pair {
 
     numeric_pair(T xp, T yp) : x(xp), y(yp) {}
 
-
     template <typename X>
     numeric_pair(const X & n) : x(n), y(0) {
     }
     
-    template <typename X>
-    numeric_pair(const numeric_pair<X> & n) : x(n.x), y(n.y) {}
+    numeric_pair(const numeric_pair<T> & n) : x(n.x), y(n.y) {}
     
     template <typename X, typename Y>
-    numeric_pair(X xp, Y yp) : numeric_pair(convert_struct<T, X>::convert(xp), convert_struct<T, Y>::convert(yp)) {}
+    numeric_pair(X xp, Y yp) : x(convert_struct<T, X>::convert(xp)), y(convert_struct<T, Y>::convert(yp)) {}
 
     bool operator<(const numeric_pair& a) const {
         return x < a.x || (x == a.x && y < a.y);
diff --git a/src/util/lp/permutation_matrix.h b/src/util/lp/permutation_matrix.h
index 6c0367482..4bdd57f25 100644
--- a/src/util/lp/permutation_matrix.h
+++ b/src/util/lp/permutation_matrix.h
@@ -132,42 +132,4 @@ class permutation_matrix : public tail_matrix<T, X> {
         
     }; // end of the permutation class
 
-#ifdef LEAN_DEBUG
-template <typename T, typename X>
-class permutation_generator {
-    unsigned m_n;
-    permutation_generator* m_lower;
-    bool m_done = false;
-    permutation_matrix<T, X> m_current;
-    unsigned m_last;
-public:
-    permutation_generator(unsigned n);
-    permutation_generator(const permutation_generator & o);
-    bool move_next();
-
-    ~permutation_generator() {
-        if (m_lower != nullptr) {
-            delete m_lower;
-        }
-    }
-
-    permutation_matrix<T, X> *current() {
-        return &m_current;
-    }
-};
-
-template <typename T, typename X>
-inline unsigned number_of_inversions(permutation_matrix<T, X> & p);
-
-template <typename T, typename X>
-int sign(permutation_matrix<T, X> & p) {
-    return is_even(number_of_inversions(p))? 1: -1;
-}
-
-template <typename T, typename X>
-T det_val_on_perm(permutation_matrix<T, X>* u, const matrix<T, X>& m);
-
-template <typename T, typename X>
-T determinant(const matrix<T, X>& m);
-#endif
 }
diff --git a/src/util/lp/permutation_matrix.hpp b/src/util/lp/permutation_matrix.hpp
index 09435674d..ec9af5a50 100644
--- a/src/util/lp/permutation_matrix.hpp
+++ b/src/util/lp/permutation_matrix.hpp
@@ -320,100 +320,4 @@ template <typename T, typename X> bool permutation_matrix<T, X>::is_identity() c
 }
 
 
-#ifdef LEAN_DEBUG
-template <typename T, typename X>
-permutation_generator<T, X>::permutation_generator(unsigned n): m_n(n), m_current(n) {
-    lean_assert(n > 0);
-    if (n > 1) {
-        m_lower = new permutation_generator(n - 1);
-    } else {
-        m_lower = nullptr;
-    }
-
-    m_last = 0;
-}
-
-template <typename T, typename X>
-permutation_generator<T, X>::permutation_generator(const permutation_generator & o): m_n(o.m_n), m_done(o.m_done), m_current(o.m_current), m_last(o.m_last) {
-    if (m_lower != nullptr) {
-        m_lower = new permutation_generator(o.m_lower);
-    } else {
-        m_lower = nullptr;
-    }
-}
-
-template <typename T, typename X> bool
-permutation_generator<T, X>::move_next() {
-    if (m_done) {
-        return false;
-    }
-
-    if (m_lower == nullptr) {
-        if (m_last == 0) {
-            m_last++;
-            return true;
-        } else {
-            m_done = true;
-            return false;
-        }
-    } else {
-        if (m_last < m_n && m_last > 0) {
-            m_current[m_last - 1] = m_current[m_last];
-            m_current[m_last] = m_n - 1;
-            m_last++;
-            return true;
-        } else {
-            if (m_lower -> move_next()) {
-                auto lower_curr = m_lower -> current();
-                for ( unsigned i = 1; i < m_n; i++ ){
-                    m_current[i] = (*lower_curr)[i - 1];
-                }
-                m_current[0] = m_n - 1;
-                m_last = 1;
-                return true;
-            } else {
-                m_done = true;
-                return false;
-            }
-        }
-    }
-}
-
-template <typename T, typename X>
-inline unsigned number_of_inversions(permutation_matrix<T, X> & p) {
-    unsigned ret = 0;
-    unsigned n = p.size();
-    for (unsigned i = 0; i < n; i++) {
-        for (unsigned j = i + 1; j < n; j++) {
-            if (p[i] > p[j]) {
-                ret++;
-            }
-        }
-    }
-    return ret;
-}
-
-template <typename T, typename X>
-T det_val_on_perm(permutation_matrix<T, X>* u, const matrix<T, X>& m) {
-    unsigned n = m.row_count();
-    T ret = numeric_traits<T>::one();
-    for (unsigned i = 0; i < n; i++) {
-        unsigned j = (*u)[i];
-        ret *=  m(i, j);
-    }
-    return ret * sign(*u);
-}
-
-template <typename T, typename X>
-T determinant(const matrix<T, X>& m) {
-    lean_assert(m.column_count() == m.row_count());
-    unsigned n = m.row_count();
-    permutation_generator<T, X> allp(n);
-    T ret = numeric_traits<T>::zero();
-    while (allp.move_next()){
-        ret += det_val_on_perm(allp.current(), m);
-    }
-    return ret;
-}
-#endif
 }
diff --git a/src/util/lp/permutation_matrix_instances.cpp b/src/util/lp/permutation_matrix_instances.cpp
index e0cc1a144..91473fabc 100644
--- a/src/util/lp/permutation_matrix_instances.cpp
+++ b/src/util/lp/permutation_matrix_instances.cpp
@@ -46,11 +46,6 @@ template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq>
 template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_reverse_from_left_to_T(vector<lean::mpq>&);
 template void lean::permutation_matrix<lean::mpq, lean::numeric_pair<lean::mpq> >::apply_reverse_from_right_to_T(vector<lean::mpq >&);
 template void lean::permutation_matrix<double, double>::multiply_by_permutation_from_right(lean::permutation_matrix<double, double>&);
-
-#ifdef LEAN_DEBUG
-template bool lean::permutation_generator<double, double>::move_next();
-template lean::permutation_generator<double, double>::permutation_generator(unsigned int);
-#endif
 template lean::permutation_matrix<double, double>::permutation_matrix(unsigned int);
 template void lean::permutation_matrix<double, double>::apply_reverse_from_left_to_X(vector<double> &);
 template void  lean::permutation_matrix< lean::mpq, lean::mpq>::apply_reverse_from_left_to_X(vector<lean::mpq> &);
diff --git a/src/util/lp/random_updater.h b/src/util/lp/random_updater.h
index 3dbab8323..2d65e191d 100644
--- a/src/util/lp/random_updater.h
+++ b/src/util/lp/random_updater.h
@@ -16,12 +16,15 @@ namespace lean {
 template <typename T> struct numeric_pair; // forward definition
 class lar_core_solver; // forward definition
 class random_updater {
-    unsigned range = 100000;
+    unsigned range ;
     struct interval {
-        bool upper_bound_is_set = false;
+        bool upper_bound_is_set;
         numeric_pair<mpq> upper_bound;
-        bool low_bound_is_set = false;
+        bool low_bound_is_set;
         numeric_pair<mpq> low_bound;
+        interval() : upper_bound_is_set(false),
+                     low_bound_is_set(false) {}
+        
         void set_low_bound(const numeric_pair<mpq> & v) {
             if (low_bound_is_set) {
                 low_bound = std::max(v, low_bound);
diff --git a/src/util/lp/random_updater.hpp b/src/util/lp/random_updater.hpp
index 68e2f5bc9..635053a09 100644
--- a/src/util/lp/random_updater.hpp
+++ b/src/util/lp/random_updater.hpp
@@ -12,7 +12,9 @@ namespace lean {
 
 random_updater::random_updater(
                                lar_core_solver & lar_core_solver,
-                               const vector<unsigned> & column_indices) : m_core_solver(lar_core_solver) {
+                               const vector<unsigned> & column_indices) :
+    m_core_solver(lar_core_solver),
+    range(100000) {
     for (unsigned j : column_indices)
         add_column_to_sets(j);
 }
diff --git a/src/util/lp/sparse_matrix.h b/src/util/lp/sparse_matrix.h
index 96f0cf2ae..8fcb75e03 100644
--- a/src/util/lp/sparse_matrix.h
+++ b/src/util/lp/sparse_matrix.h
@@ -30,10 +30,10 @@ class sparse_matrix
 #endif
 {
     struct col_header {
-        unsigned m_shortened_markovitz = 0;
+        unsigned m_shortened_markovitz;
         vector<indexed_value<T>> m_values; // the actual column values
 
-        col_header()  {}
+        col_header(): m_shortened_markovitz(0) {}
 
         void shorten_markovich_by_one() {
             m_shortened_markovitz++;
@@ -44,7 +44,7 @@ class sparse_matrix
         }
     };
 
-    unsigned m_n_of_active_elems = 0;
+    unsigned m_n_of_active_elems;
     binary_heap_upair_queue<unsigned> m_pivot_queue;
 public:
     vector<vector<indexed_value<T>>>  m_rows;
diff --git a/src/util/lp/sparse_matrix.hpp b/src/util/lp/sparse_matrix.hpp
index 0d2a90ea0..f25eadbb4 100644
--- a/src/util/lp/sparse_matrix.hpp
+++ b/src/util/lp/sparse_matrix.hpp
@@ -40,7 +40,8 @@ sparse_matrix<T, X>::sparse_matrix(static_matrix<T, X> const &A, vector<unsigned
     m_row_permutation(A.row_count()),
     m_column_permutation(A.row_count()),
     m_work_pivot_vector(A.row_count(), -1),
-    m_processed(A.row_count()) {
+    m_processed(A.row_count()),
+    m_n_of_active_elems(0) {
     init_row_headers();
     init_column_headers();
     copy_B(A, basis);
diff --git a/src/util/lp/sparse_vector.h b/src/util/lp/sparse_vector.h
index 78d7ff5be..975cb7f28 100644
--- a/src/util/lp/sparse_vector.h
+++ b/src/util/lp/sparse_vector.h
@@ -20,7 +20,7 @@ public:
     }
 #ifdef LEAN_DEBUG
     T operator[] (unsigned i) const {
-        for (auto t : m_data) {
+        for (auto &t : m_data) {
             if (t.first == i) return t.second;
         }
         return numeric_traits<T>::zero();
diff --git a/src/util/lp/square_dense_submatrix.h b/src/util/lp/square_dense_submatrix.h
index 10ae973d6..019497aa5 100644
--- a/src/util/lp/square_dense_submatrix.h
+++ b/src/util/lp/square_dense_submatrix.h
@@ -42,7 +42,7 @@ public:
     unsigned m_index_start;
     unsigned m_dim;
     vector<T> m_v;
-    sparse_matrix<T, X> * m_parent = nullptr;
+    sparse_matrix<T, X> * m_parent;
     permutation_matrix<T, X>  m_row_permutation;
     indexed_vector<T> m_work_vector;
 public:
diff --git a/src/util/lp/ul_pair.h b/src/util/lp/ul_pair.h
index 0e96364ce..2e77a7db0 100644
--- a/src/util/lp/ul_pair.h
+++ b/src/util/lp/ul_pair.h
@@ -32,14 +32,14 @@ inline bool compare(const std::pair<mpq, var_index> & a, const std::pair<mpq, va
 }
 
 class ul_pair {
-    constraint_index m_low_bound_witness = static_cast<constraint_index>(-1);
-    constraint_index m_upper_bound_witness = static_cast<constraint_index>(-1);
+    constraint_index m_low_bound_witness;
+    constraint_index m_upper_bound_witness;
 public:
     constraint_index& low_bound_witness() {return m_low_bound_witness;}
     constraint_index low_bound_witness() const {return m_low_bound_witness;}
     constraint_index& upper_bound_witness() { return m_upper_bound_witness;}
     constraint_index upper_bound_witness() const {return m_upper_bound_witness;}
-    row_index m_i = static_cast<row_index>(-1);
+    row_index m_i;
     bool operator!=(const ul_pair & p) const {
         return !(*this == p);
     }
@@ -50,8 +50,15 @@ public:
             m_i == p.m_i;
     }
 	// empty constructor
-	ul_pair(){}
-    ul_pair(row_index ri) : m_i(ri) {}
+	ul_pair() :
+        m_low_bound_witness(static_cast<constraint_index>(-1)),
+        m_upper_bound_witness(static_cast<constraint_index>(-1)),
+        m_i(static_cast<row_index>(-1))
+{}
+    ul_pair(row_index ri) :
+        m_low_bound_witness(static_cast<constraint_index>(-1)),
+        m_upper_bound_witness(static_cast<constraint_index>(-1)),
+        m_i(ri) {}
     ul_pair(const ul_pair & o): m_low_bound_witness(o.m_low_bound_witness), m_upper_bound_witness(o.m_upper_bound_witness), m_i(o.m_i) {}
 };