diff --git a/src/util/lp/lp_core_solver_base.cpp b/src/util/lp/lp_core_solver_base.cpp
index d8cf502856..d995aa3007 100644
--- a/src/util/lp/lp_core_solver_base.cpp
+++ b/src/util/lp/lp_core_solver_base.cpp
@@ -198,11 +198,11 @@ A_mult_x_is_off() {
         X eps = feps * (one + T(0.1) * abs(m_b[i]));
 
         if (delta > eps) {
-            std::cout << "x is off (";
-            std::cout << "m_b[" << i  << "] = " << m_b[i] << " ";
-            std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
-            std::cout << "delta = " << delta << ' ';
-            std::cout << "iters = " << m_total_iterations << ")" << std::endl;
+            // std::cout << "x is off (";
+            // std::cout << "m_b[" << i  << "] = " << m_b[i] << " ";
+            // std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
+            // std::cout << "delta = " << delta << ' ';
+            // std::cout << "iters = " << m_total_iterations << ")" << std::endl;
             return true;
         }
     }
@@ -333,7 +333,7 @@ set_non_basic_x_to_correct_bounds() {
     }
 }
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
-column_is_dual_feasible(unsigned j) {
+column_is_dual_feasible(unsigned j) const {
     switch (m_column_type[j]) {
     case fixed:
     case boxed:
@@ -353,7 +353,7 @@ column_is_dual_feasible(unsigned j) {
     }
 }
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
-d_is_not_negative(unsigned j) {
+d_is_not_negative(unsigned j) const {
     if (numeric_traits<T>::precise()) {
         return m_d[j] >= numeric_traits<T>::zero();
     }
@@ -361,7 +361,7 @@ d_is_not_negative(unsigned j) {
 }
 
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
-d_is_not_positive(unsigned j) {
+d_is_not_positive(unsigned j) const {
     if (numeric_traits<T>::precise()) {
         return m_d[j] <= numeric_traits<T>::zero();
     }
@@ -686,4 +686,29 @@ get_non_basic_column_value_position(unsigned j) {
         lean_unreachable();
     }
 }
+
+template <typename T, typename X> void lp_core_solver_base<T, X>::init_lu() {
+    init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
+    this->m_refactor_counter = 0;
 }
+
+template <typename T, typename X> int lp_core_solver_base<T, X>::pivots_in_column_and_row_are_different(int entering, int leaving) const {
+    const T & column_p = this->m_ed[this->m_basis_heading[leaving]];
+    const T & row_p = this->m_pivot_row[entering];
+    if (is_zero(column_p) || is_zero(row_p)) return true; // pivots cannot be zero
+    // the pivots have to have the same sign
+    if (column_p < 0) {
+        if (row_p > 0)
+            return 2;
+    } else { // column_p > 0
+        if (row_p < 0)
+            return 2;
+    }
+    T diff_normalized = abs((column_p - row_p) / (numeric_traits<T>::one() + abs(row_p)));
+    if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10)))
+        return 1;
+    return 0;
+}
+
+}
+
diff --git a/src/util/lp/lp_core_solver_base.h b/src/util/lp/lp_core_solver_base.h
index 3b8102d962..cc0c634d63 100644
--- a/src/util/lp/lp_core_solver_base.h
+++ b/src/util/lp/lp_core_solver_base.h
@@ -163,6 +163,15 @@ public:
         return below_bound(m_x[p], m_low_bound_values[p]);
     }
 
+    bool x_above_low_bound(unsigned p) {
+        return above_bound(m_x[p], m_low_bound_values[p]);
+    }
+
+    bool x_below_upper_bound(unsigned p) {
+        return below_bound(m_x[p], m_upper_bound_values[p]);
+    }
+
+    
     bool x_above_upper_bound(unsigned p) {
         return above_bound(m_x[p], m_upper_bound_values[p]);
     }
@@ -177,10 +186,11 @@ public:
         return x_is_at_low_bound(j) || x_is_at_upper_bound(j);
     }
 
-    bool column_is_dual_feasible(unsigned j);
-    bool d_is_not_negative(unsigned j);
+    bool column_is_dual_feasible(unsigned j) const;
 
-    bool d_is_not_positive(unsigned j);
+    bool d_is_not_negative(unsigned j) const;
+
+    bool d_is_not_positive(unsigned j) const;
 
 
     bool time_is_over();
@@ -191,7 +201,6 @@ public:
 
     bool update_basis_and_x(int entering, int leaving, X const & tt);
 
-
     void init_basis_heading();
 
     bool basis_has_no_doubles();
@@ -244,5 +253,8 @@ public:
     void init_reduced_costs_for_one_iteration();
 
     non_basic_column_value_position get_non_basic_column_value_position(unsigned j);
+
+    void init_lu();
+    int pivots_in_column_and_row_are_different(int entering, int leaving) const;
 };
 }
diff --git a/src/util/lp/lp_core_solver_base_instances.cpp b/src/util/lp/lp_core_solver_base_instances.cpp
index 3c54e6df13..a3ec83e663 100644
--- a/src/util/lp/lp_core_solver_base_instances.cpp
+++ b/src/util/lp/lp_core_solver_base_instances.cpp
@@ -9,7 +9,7 @@ template bool lean::lp_core_solver_base<double, double>::A_mult_x_is_off();
 template bool lean::lp_core_solver_base<double, double>::basis_heading_is_correct();
 template void lean::lp_core_solver_base<double, double>::calculate_pivot_row_of_B_1(unsigned int);
 template void lean::lp_core_solver_base<double, double>::calculate_pivot_row_when_pivot_row_of_B1_is_ready();
-template bool lean::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int);
+template bool lean::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int) const;
 template void lean::lp_core_solver_base<double, double>::fill_reduced_costs_from_m_y_by_rows();
 template bool lean::lp_core_solver_base<double, double>::find_x_by_solving();
 template lean::non_basic_column_value_position lean::lp_core_solver_base<double, double>::get_non_basic_column_value_position(unsigned int);
@@ -29,7 +29,7 @@ template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::A_mult_x_is_off()
 template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::basis_heading_is_correct();
 template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::calculate_pivot_row_of_B_1(unsigned int);
 template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::calculate_pivot_row_when_pivot_row_of_B1_is_ready();
-template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::column_is_dual_feasible(unsigned int);
+template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::column_is_dual_feasible(unsigned int) const;
 template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::fill_reduced_costs_from_m_y_by_rows();
 template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::find_x_by_solving();
 template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::init_reduced_costs_for_one_iteration();
@@ -67,3 +67,7 @@ template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq>
 template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::restore_state(lean::mpq*, lean::mpq*);
 template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::save_state(lean::mpq*, lean::mpq*);
 template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::solve_yB(std::vector<lean::mpq, std::allocator<lean::mpq> >&);
+template void lean::lp_core_solver_base<double, double>::init_lu();
+template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::init_lu();
+template int lean::lp_core_solver_base<double, double>::pivots_in_column_and_row_are_different(int, int) const;
+template int lean::lp_core_solver_base<lean::mpq, lean::mpq>::pivots_in_column_and_row_are_different(int, int) const;
diff --git a/src/util/lp/lp_dual_core_solver.cpp b/src/util/lp/lp_dual_core_solver.cpp
index d34a87a3af..f6c796f922 100644
--- a/src/util/lp/lp_dual_core_solver.cpp
+++ b/src/util/lp/lp_dual_core_solver.cpp
@@ -80,29 +80,22 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_non_ba
 
 template <typename T, typename X> bool lp_dual_core_solver<T, X>::update_basis(int entering, int leaving) {
     // the second argument is the element of the entering column from the pivot row - its value should be equal to the low diagonal element of the bump after all pivoting is done
-    if (!(this->m_refactor_counter++ >= 200)) {
+    if (this->m_refactor_counter++ < 200) {
         this->m_factorization->replace_column(leaving, this->m_ed[this->m_factorization->basis_heading(leaving)], this->m_w);
-        if (this->m_factorization->get_status() != LU_status::OK) {
-#ifdef LEAN_DEBUG
-            std::cout << "failed on replace_column( " << leaving << ", " <<  this->m_ed[this->m_factorization->basis_heading(leaving)] << ") total_iterations = " << this->m_total_iterations << std::endl;
-#endif
-            init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
-            this->m_iters_with_no_cost_growing++;
-            this->m_status = UNSTABLE;
-            return false;
-        }
-        this->m_factorization->change_basis(entering, leaving);
-    } else { // need to refactor
-        this->m_factorization->change_basis(entering, leaving);
-        init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
-        if (this->m_factorization->get_status() != LU_status::OK) {
-#ifdef LEAN_DEBUG
-            std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << std::endl;
-#endif
-            this->m_iters_with_no_cost_growing++;
-            return false;
+        if (this->m_factorization->get_status() == LU_status::OK) {
+            this->m_factorization->change_basis(entering, leaving);
+            return true;
         }
     }
+    // need to refactor
+    this->m_factorization->change_basis(entering, leaving);
+    init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
+    this->m_refactor_counter = 0;
+    if (this->m_factorization->get_status() != LU_status::OK) {
+        std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << std::endl;
+        this->m_iters_with_no_cost_growing++;
+        return false;
+    }
     return true;
 }
 
@@ -204,19 +197,24 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
     unsigned i = offset_in_rows;
     unsigned rows_left = number_of_rows_to_try;
     do  {
-    loop_start:
+        if (m_forbidden_rows.find(i) != m_forbidden_rows.end()) {
+            if (++i == this->m_m) {
+                i = 0;
+            }
+            continue;
+        }
         T se = pricing_for_row(i);
         if (se > steepest_edge_max) {
             steepest_edge_max = se;
             m_r = i;
+            if (rows_left > 0) {
+                rows_left--;
+            }
         }
         if (++i == this->m_m) {
             i = 0;
         }
-        if (rows_left > 0) {
-            rows_left--;
-        }
-    }  while (rows_left || (steepest_edge_max < T(1e-5) && i != initial_offset_in_rows));
+    } while (i != initial_offset_in_rows && (rows_left || steepest_edge_max < T(1e-5)));
     if (m_r == -1) {
         if (this->m_status != UNSTABLE) {
             this->m_status = OPTIMAL;
@@ -225,22 +223,24 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
         m_p = this->m_basis[m_r];
         m_delta = get_delta();
         if (advance_on_known_p()){
+            m_forbidden_rows.clear();
             return;
         }
+        // failure in advance_on_known_p
         if (this->m_status == FLOATING_POINT_ERROR) {
             return;
         }
         this->m_status = UNSTABLE;
-        if (i == initial_offset_in_rows) {
-            return; // made a full loop
-        }
-        m_p = m_r = -1;
-        steepest_edge_max = numeric_traits<T>::zero();
-        rows_left = number_of_rows_to_try;
-        goto loop_start;
+        m_forbidden_rows.insert(m_r);
     }
 }
 
+    // this calculation is needed for the steepest edge update,
+    // it hijackes m_pivot_row_of_B_1 for this purpose since we will need it anymore to the end of the cycle
+template <typename T, typename X> void lp_dual_core_solver<T, X>::DSE_FTran() { // todo, see algorithm 7 from page 35
+    this->m_factorization->solve_By(this->m_pivot_row_of_B_1);
+}
+
 template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_known_p() {
     if (done()) {
         return true;
@@ -251,7 +251,15 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_kno
         return true;
     }
     calculate_beta_r_precisely();
-    FTran();
+    this->solve_Bd(m_q); //FTRAN
+    int pivot_compare_result = this->pivots_in_column_and_row_are_different(m_q, m_p);
+    if (!pivot_compare_result){;}
+    else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
+        lean_unreachable(); // not implemented yet
+    } else {
+        lean_assert(pivot_compare_result == 1);
+        this->init_lu();
+    }
     DSE_FTran();
     return basis_change_and_update();
 }
@@ -315,9 +323,9 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::fill_breakpoin
     }
 }
 
-template <typename T, typename X> void lp_dual_core_solver<T, X>::FTran() {
-    this->solve_Bd(m_q);
-}
+// template <typename T, typename X> void lp_dual_core_solver<T, X>::FTran() {
+//     this->solve_Bd(m_q);
+// }
 
 template <typename T, typename X> T lp_dual_core_solver<T, X>::get_delta() {
     switch (this->m_column_type[m_p]) {
@@ -451,46 +459,115 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::init_betas_pre
 template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_and_update() {
     update_betas();
     update_d_and_xB();
-    m_theta_P = m_delta / this->m_ed[m_r];
-    xb_minus_delta_p_pivot_column();
+    //    m_theta_P = m_delta / this->m_ed[m_r];
+    m_theta_P = m_delta / this->m_pivot_row[m_q];
+    //    xb_minus_delta_p_pivot_column();
     apply_flips();
-    this->m_x[m_q] += m_theta_P;
-    if (!update_basis(m_q, m_p) || this->A_mult_x_is_off() || !problem_is_dual_feasible()) {
+    if (!this->update_basis_and_x(m_q, m_p, m_theta_P)) {
+          init_betas_precisely();
+          return false;
+    }
+
+    if (snap_runaway_nonbasic_column(m_p)) {
+        if(!this->find_x_by_solving()) {
+            revert_to_previous_basis();
+            this->m_iters_with_no_cost_growing++;
+            return false;
+        }
+    }
+
+    if (!problem_is_dual_feasible() ) {
+        // todo : shift the costs!!!!
         revert_to_previous_basis();
         this->m_iters_with_no_cost_growing++;
         return false;
     }
+    
     lean_assert(d_is_correct());
     return true;
 }
 
+template <typename T, typename X> void lp_dual_core_solver<T, X>::recover_leaving() {
+    switch (m_entering_boundary_position) {
+    case at_low_bound:
+    case at_fixed:
+        this->m_x[m_q] = this->m_low_bound_values[m_q];
+        break;
+    case at_upper_bound:
+        this->m_x[m_q] = this->m_upper_bound_values[m_q];
+        break;
+    case free_of_bounds:
+        this->m_x[m_q] = zero_of_type<X>();
+    default:
+        lean_unreachable();
+    }
+}
+
 template <typename T, typename X> void lp_dual_core_solver<T, X>::revert_to_previous_basis() {
+    //    std::cout << "recovering basis p = " << m_p << " q = " << m_q << std::endl;
     this->m_factorization->change_basis(m_p, m_q);
     init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
     if (this->m_factorization->get_status() != LU_status::OK) {
+        this->m_status = FLOATING_POINT_ERROR; // complete failure
+        return;
+    }
+    recover_leaving();
+    if (!this->find_x_by_solving()) {
         this->m_status = FLOATING_POINT_ERROR;
         return;
     }
-    this->m_x[m_q] -= m_theta_P;
-    snap_xN_to_bounds();
     recalculate_xB_and_d();
-    if (this->A_mult_x_is_off()) {
-        this->m_status = FLOATING_POINT_ERROR;
-        return;
-    }
     init_betas_precisely();
 }
 
+// returns true if the column has been snapped
+template <typename T, typename X> bool lp_dual_core_solver<T, X>::snap_runaway_nonbasic_column(unsigned j) {
+    switch (this->m_column_type[j]) {
+    case fixed:
+    case low_bound:
+        if (! this->x_is_at_low_bound(j)) { 
+            this->m_x[j] = this->m_low_bound_values[j];
+            return true;
+        }
+        break;
+    case boxed:
+        {
+            bool closer_to_low_bound = abs(this->m_low_bound_values[j] - this->m_x[j]) < abs(this->m_upper_bound_values[j] - this->m_x[j]);
+            if (closer_to_low_bound) {
+                if (! this->x_is_at_low_bound(j)) { 
+                    this->m_x[j] = this->m_low_bound_values[j];
+                    return true;
+                }
+            } else {
+                if (! this->x_is_at_upper_bound(j)) { 
+                    this->m_x[j] = this->m_low_bound_values[j];
+                    return true;
+                }
+            }
+        }
+        break;
+    case upper_bound:
+        if (! this->x_is_at_upper_bound(j)) { 
+            this->m_x[j] = this->m_upper_bound_values[j];
+            return true;
+        }
+        break;
+    default:
+        break;
+    }
+    return false;
+}
 
-template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() {
+
+template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() const {
     for (unsigned j : non_basis()){
         if (!this->column_is_dual_feasible(j)) {
-            std::cout << "column " << j << " is not dual feasible" << std::endl;
-            std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
-            std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
-            std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
-            std::cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << std::endl;
-            std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
+            // std::cout << "column " << j << " is not dual feasible" << std::endl;
+            // std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
+            // std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
+            // std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
+            // std::cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << std::endl;
+            // std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
             return false;
         }
     }
@@ -500,7 +577,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dua
 template <typename T, typename X> unsigned lp_dual_core_solver<T, X>::get_number_of_rows_to_try_for_leaving() {
     unsigned s = this->m_m;
     if (this->m_m > 300) {
-        s = (unsigned)(s / this->m_settings.percent_of_entering_to_check * 100);
+        s = (unsigned)((s / 100.0) * this->m_settings.percent_of_entering_to_check);
     }
     return my_random() % s + 1;
 }
diff --git a/src/util/lp/lp_dual_core_solver.h b/src/util/lp/lp_dual_core_solver.h
index ed05226da1..87bbcfac01 100644
--- a/src/util/lp/lp_dual_core_solver.h
+++ b/src/util/lp/lp_dual_core_solver.h
@@ -32,7 +32,7 @@ public:
     std::vector<T> m_a_wave;
     std::vector<T> m_betas; // m_betas[i] is approximately a square of the norm of the i-th row of the reverse of B
     T m_harris_tolerance;
-
+    std::set<unsigned> m_forbidden_rows;
     lp_dual_core_solver(static_matrix<T, X> & A,
                         std::vector<bool> & can_enter_basis,
                         std::vector<X> & b, // the right side std::vector
@@ -59,7 +59,7 @@ public:
 
     void recalculate_d();
 
-    std::vector<unsigned> & non_basis() { return this->m_factorization->m_non_basic_columns; }
+    const std::vector<unsigned> & non_basis() const { return this->m_factorization->m_non_basic_columns; }
 
     void init_betas();
 
@@ -85,14 +85,7 @@ public:
 
     void fill_breakpoint_set();
 
-    void FTran();
-
-    // this calculation is needed for the steepest edge update,
-    // it hijackes m_pivot_row_of_B_1 for this purpose since we will need it anymore to the end of the cycle
-    void DSE_FTran() { // todo, see algorithm 7 from page 35
-        this->m_factorization->solve_By(this->m_pivot_row_of_B_1);
-    }
-
+    void DSE_FTran();
     T get_delta();
 
     void restore_d();
@@ -118,8 +111,16 @@ public:
 
     void revert_to_previous_basis();
 
-    bool problem_is_dual_feasible();
+    non_basic_column_value_position m_entering_boundary_position;
+    bool update_basis_and_x_local(int entering, int leaving, X const & tt);
+    void recover_leaving();
 
+    bool problem_is_dual_feasible() const;
+
+    bool snap_runaway_nonbasic_column(unsigned);
+
+    bool snap_runaway_nonbasic_columns();
+    
     unsigned get_number_of_rows_to_try_for_leaving();
 
     void update_a_wave(const T & del, unsigned j) {
diff --git a/src/util/lp/lp_dual_core_solver_instances.cpp b/src/util/lp/lp_dual_core_solver_instances.cpp
index 46306fed00..378b21962f 100644
--- a/src/util/lp/lp_dual_core_solver_instances.cpp
+++ b/src/util/lp/lp_dual_core_solver_instances.cpp
@@ -13,3 +13,5 @@ template void lean::lp_dual_core_solver<double, double>::start_with_initial_basi
 template void lean::lp_dual_core_solver<double, double>::solve();
 template void lean::lp_dual_core_solver<lean::mpq, lean::mpq>::restore_non_basis();
 template void lean::lp_dual_core_solver<double, double>::restore_non_basis();
+template void lean::lp_dual_core_solver<double, double>::revert_to_previous_basis();
+template void lean::lp_dual_core_solver<lean::mpq, lean::mpq>::revert_to_previous_basis();
diff --git a/src/util/lp/lp_dual_simplex.cpp b/src/util/lp/lp_dual_simplex.cpp
index 846803d965..c16a59a2cb 100644
--- a/src/util/lp/lp_dual_simplex.cpp
+++ b/src/util/lp/lp_dual_simplex.cpp
@@ -120,6 +120,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2()
         lean_unreachable();
     }
     this->m_second_stage_iterations = m_core_solver->m_total_iterations;
+    this->m_total_iterations = this->m_first_stage_iterations + this->m_second_stage_iterations;
 }
 
 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros() {
@@ -272,7 +273,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_s
             this->m_costs[artificial] = numeric_traits<T>::zero();
             artificial++;
         } else {
-            // we can put a slack_var into the basis, and atemplate <typename T, typename X> void lp_dual_simplex<T, X>::adding an artificial variable
+            // we can put a slack_var into the basis, and avoid adding an artificial variable
             this->m_basis[row] = slack_var;
             this->m_costs[slack_var] = numeric_traits<T>::zero();
         }
@@ -290,7 +291,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_s
             this->m_costs[artificial] = numeric_traits<T>::zero();
             artificial++;
         } else {
-            // we can put slack_var into the basis, and atemplate <typename T, typename X> void lp_dual_simplex<T, X>::adding an artificial variable
+            // we can put slack_var into the basis, and avoid adding an artificial variable
             this->m_basis[row] = slack_var;
             this->m_costs[slack_var] = numeric_traits<T>::zero();
         }
diff --git a/src/util/lp/lp_primal_core_solver.cpp b/src/util/lp/lp_primal_core_solver.cpp
index ff5b618668..3dfda91f31 100644
--- a/src/util/lp/lp_primal_core_solver.cpp
+++ b/src/util/lp/lp_primal_core_solver.cpp
@@ -261,11 +261,6 @@ lp_primal_core_solver(static_matrix<T, X> & A,
 #endif
 }
 
-template <typename T, typename X>    void lp_primal_core_solver<T, X>::init_lu() {
-    init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
-    this->m_refactor_counter = 0;
-}
-
 template <typename T, typename X> bool lp_primal_core_solver<T, X>::initial_x_is_correct() {
     std::set<unsigned> basis_set;
     for (int i = 0; i < this->m_A.row_count(); i++) {
@@ -396,23 +391,6 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::init_run(
     init_reduced_costs();
 }
 
-template <typename T, typename X>    int lp_primal_core_solver<T, X>::pivots_in_column_and_row_are_different(int entering, int leaving) const {
-    const T & column_p = this->m_ed[this->m_basis_heading[leaving]];
-    const T & row_p = this->m_pivot_row[entering];
-    if (is_zero(column_p) || is_zero(row_p)) return true; // pivots cannot be zero
-    // the pivots have to have the same sign
-    if (column_p < 0) {
-        if (row_p > 0)
-            return 2;
-    } else { // column_p > 0
-        if (row_p < 0)
-            return 2;
-    }
-    T diff_normalized = abs((column_p - row_p) / (numeric_traits<T>::one() + abs(row_p)));
-    if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10)))
-        return 1;
-    return 0;
-}
 
 template <typename T, typename X>    void lp_primal_core_solver<T, X>::calc_working_vector_beta_for_column_norms(){
     unsigned i = this->m_m;
@@ -429,7 +407,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
         lean_assert(!this->A_mult_x_is_off() );
         this->update_x(entering, t * m_sign_of_entering_delta);
         if (this->A_mult_x_is_off() && !this->find_x_by_solving()) {
-            init_lu();
+            this->init_lu();
             if (!this->find_x_by_solving()) {
                 this->restore_x(entering, t * m_sign_of_entering_delta);
                 m_forbidden_enterings.insert(entering);
@@ -447,7 +425,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
     unsigned pivot_row = this->m_factorization->basis_heading(leaving);
     this->calculate_pivot_row_of_B_1(pivot_row);
     this->calculate_pivot_row_when_pivot_row_of_B1_is_ready();
-    int pivot_compare_result = pivots_in_column_and_row_are_different(entering, leaving);
+    int pivot_compare_result = this->pivots_in_column_and_row_are_different(entering, leaving);
     if (!pivot_compare_result){;}
     else if (pivot_compare_result == 2) { // the sign is changed, cannot continue
         m_forbidden_enterings.insert(entering);
@@ -456,7 +434,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
         return;
     } else {
         lean_assert(pivot_compare_result == 1);
-        init_lu();
+        this->init_lu();
     }
     calc_working_vector_beta_for_column_norms();
     if (!this->update_basis_and_x(entering, leaving, t * m_sign_of_entering_delta)) {
@@ -491,7 +469,7 @@ template <typename T, typename X>    void lp_primal_core_solver<T, X>::advance_o
     this->solve_Bd(entering);
     int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
     if (refresh_result) {
-        init_lu();
+        this->init_lu();
         init_reduced_costs();
         if (refresh_result == 2) {
             this->m_iters_with_no_cost_growing++;
@@ -565,7 +543,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
         case OPTIMAL:  // double check that we are at optimum
         case INFEASIBLE:
             m_forbidden_enterings.clear();
-            init_lu();
+            this->init_lu();
             lean_assert(this->m_factorization->get_status() == LU_status::OK);
             set_current_x_is_feasible();
             init_reduced_costs();
@@ -577,7 +555,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
             break;
         case TENTATIVE_UNBOUNDED:
             m_forbidden_enterings.clear();
-            init_lu();
+            this->init_lu();
             lean_assert(this->m_factorization->get_status() == LU_status::OK);
             init_reduced_costs();
             break;
@@ -586,7 +564,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
 
         case UNSTABLE:
             // m_forbidden_enterings.clear();
-            init_lu();
+            this->init_lu();
             init_reduced_costs();
             this->m_status = UNKNOWN;
             break;
diff --git a/src/util/lp/lp_primal_core_solver.h b/src/util/lp/lp_primal_core_solver.h
index 221d00a5be..4629e54d1b 100644
--- a/src/util/lp/lp_primal_core_solver.h
+++ b/src/util/lp/lp_primal_core_solver.h
@@ -127,7 +127,6 @@ public:
                           std::vector<X> & upper_bound_values,
                           lp_settings & settings,
                           std::unordered_map<unsigned, std::string> const & column_names);
-    void init_lu();
     bool initial_x_is_correct();
 #ifdef LEAN_DEBUG
     void check_Ax_equal_b();
@@ -148,8 +147,6 @@ public:
 
     void init_run();
 
-    int pivots_in_column_and_row_are_different(int entering, int leaving) const;
-
     void calc_working_vector_beta_for_column_norms();
 
     void advance_on_entering_and_leaving(int entering, int leaving, X & t);