remove m_b from lar_core_solver

the column vector is pure overhead for the way the lar solver uses lp.
Some other solver modules use column vectors b and integrate with the lp_core_solver_base. The interaction model should be reviewed.
Unused solvers should be removed to make it easier to maintain this code.

src/math/lp/lar_core_solver.h

@@ -66,6 +66,7 @@ public:
                     );
 
     lp_settings & settings() { return m_r_solver.m_settings;}
+
     const lp_settings & settings() const { return m_r_solver.m_settings;}
     
     int get_infeasible_sum_sign() const { return m_infeasible_sum_sign;   }
@@ -340,6 +341,7 @@ public:
                   switch (m_column_types[j]) {
                   case column_type::free_column:
                       lp_assert(false); // unreachable
+                      break;
                   case column_type::upper_bound:
                       s.m_x[j] = s.m_upper_bounds[j];
                       break;
@@ -365,7 +367,7 @@ public:
             }
         }
 
-        lp_assert(is_zero_vector(s.m_b));
+        // lp_assert(is_zero_vector(s.m_b));
         s.solve_Ax_eq_b();
     }
 
@@ -463,7 +465,8 @@ public:
             m_d_nbasis = m_r_nbasis;
             delete m_d_solver.m_factorization;
             m_d_solver.m_factorization = nullptr;
-        } else {
+        } 
+        else {
             prepare_solver_x_with_signature_tableau(signature);
             m_r_solver.start_tracing_basis_changes();
             m_r_solver.find_feasible_solution();

src/math/lp/lar_core_solver_def.h

@@ -46,14 +46,10 @@ lar_core_solver::lar_core_solver(
                column_names) {
 }    
     
-    
-
 void lar_core_solver::calculate_pivot_row(unsigned i) {
     m_r_solver.calculate_pivot_row(i);
 }
 
-
-
 void lar_core_solver::prefix_r() {
     if (!m_r_solver.m_settings.use_tableau()) {
         m_r_solver.m_copy_of_xB.resize(m_r_solver.m_n());
@@ -67,7 +63,7 @@ void lar_core_solver::prefix_r() {
         init_column_row_nz_for_r_solver();
     }
 
-    m_r_solver.m_b.resize(m_r_solver.m_m());
+    // m_r_solver.m_b.resize(m_r_solver.m_m());
     if (m_r_solver.m_settings.simplex_strategy() != simplex_strategy_enum::tableau_rows) {
         if(m_r_solver.m_settings.use_breakpoints_in_feasibility_search)
             m_r_solver.m_breakpoint_indices_queue.resize(m_r_solver.m_n());
@@ -78,7 +74,7 @@ void lar_core_solver::prefix_r() {
 }
 
 void lar_core_solver::prefix_d() {
-    m_d_solver.m_b.resize(m_d_solver.m_m());
+    // m_d_solver.m_b.resize(m_d_solver.m_m());
     m_d_solver.m_breakpoint_indices_queue.resize(m_d_solver.m_n());
     m_d_solver.m_copy_of_xB.resize(m_d_solver.m_n());
     m_d_solver.m_costs.resize(m_d_solver.m_n());
@@ -100,9 +96,8 @@ void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() {
     unsigned bj = m_r_basis[m_r_solver.m_inf_row_index_for_tableau];
     m_infeasible_sum_sign =  m_r_solver.inf_sign_of_column(bj);
     m_infeasible_linear_combination.clear();
-    for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) {
+    for (auto & rc : m_r_solver.m_A.m_rows[m_r_solver.m_inf_row_index_for_tableau]) 
         m_infeasible_linear_combination.push_back(std::make_pair(rc.coeff(), rc.var()));    
-    }
 }
 
 void lar_core_solver::fill_not_improvable_zero_sum() {
@@ -115,26 +110,23 @@ void lar_core_solver::fill_not_improvable_zero_sum() {
     m_infeasible_linear_combination.clear();
     for (auto j : m_r_solver.m_basis) {
         const mpq & cost_j = m_r_solver.m_costs[j];
-        if (!numeric_traits<mpq>::is_zero(cost_j)) {
+        if (!numeric_traits<mpq>::is_zero(cost_j)) 
             m_infeasible_linear_combination.push_back(std::make_pair(cost_j, j));        
-        }
     }
     // m_costs are expressed by m_d ( additional costs), substructing the latter gives 0
     for (unsigned j = 0; j < m_r_solver.m_n(); j++) {
         if (m_r_solver.m_basis_heading[j] >= 0) continue;
         const mpq & d_j = m_r_solver.m_d[j];
-        if (!numeric_traits<mpq>::is_zero(d_j)) {
+        if (!numeric_traits<mpq>::is_zero(d_j)) 
             m_infeasible_linear_combination.push_back(std::make_pair(-d_j, j));        
-        }
     }
 }
 
 unsigned lar_core_solver::get_number_of_non_ints() const {
     unsigned n = 0;
-    for (auto & x : m_r_solver.m_x) {
+    for (auto & x : m_r_solver.m_x) 
-        if (x.is_int() == false)
+        if (!x.is_int())
             n++;    
-    }
     return n;
 }
 

src/math/lp/lar_solver.cpp

@@ -8,8 +8,8 @@
 
 namespace lp {
 
-    ////////////////// methods ////////////////////////////////
     static_matrix<double, double>& lar_solver::A_d() { return m_mpq_lar_core_solver.m_d_A; }
+
     static_matrix<double, double > const& lar_solver::A_d() const { return m_mpq_lar_core_solver.m_d_A; }
 
     lp_settings& lar_solver::settings() { return m_settings; }
@@ -18,7 +18,6 @@ namespace lp {
 
     statistics& lar_solver::stats() { return m_settings.stats(); }
 
-
     void lar_solver::updt_params(params_ref const& _p) {
         smt_params_helper p(_p);
         set_track_pivoted_rows(p.arith_bprop_on_pivoted_rows());
@@ -42,7 +41,6 @@ namespace lp {
     }
 
     lar_solver::~lar_solver() {
-
         for (auto t : m_terms)
             delete t;
     }
@@ -1406,7 +1404,6 @@ namespace lp {
         A_r().m_rows.pop_back();
         A_r().m_columns.pop_back();
         CASSERT("check_static_matrix", A_r().is_correct());
-        slv.m_b.pop_back();
     }
 
     void lar_solver::remove_last_column_from_A() {
@@ -1716,8 +1713,6 @@ namespace lp {
         m_terms.push_back(t);
     }
 
-
-
     // terms
     bool lar_solver::all_vars_are_registered(const vector<std::pair<mpq, var_index>>& coeffs) {
         for (const auto& p : coeffs) {
@@ -1788,7 +1783,6 @@ namespace lp {
             A_r().fill_last_row_with_pivoting(*term,
                 j,
                 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);

src/math/lp/lp_core_solver_base.cpp

@@ -36,8 +36,8 @@ template lp::non_basic_column_value_position lp::lp_core_solver_base<lp::mpq, lp
 template lp::non_basic_column_value_position lp::lp_core_solver_base<lp::mpq, lp::mpq>::get_non_basic_column_value_position(unsigned int) const;
 template void lp::lp_core_solver_base<double, double>::init_reduced_costs_for_one_iteration();
 template lp::lp_core_solver_base<double, double>::lp_core_solver_base(
-    lp::static_matrix<double, double>&, vector<double>&, 
+    lp::static_matrix<double, double>&, // vector<double>&, 
-    vector<unsigned int >&,
+    vector<unsigned>&,
     vector<unsigned> &, vector<int> &,
     vector<double >&, 
     vector<double >&, 
@@ -80,7 +80,9 @@ template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calc
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init();
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init_basis_heading_and_non_basic_columns_vector();
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init_reduced_costs_for_one_iteration();
-template lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::lp_core_solver_base(lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, vector<lp::numeric_pair<lp::mpq> >&, vector<unsigned int >&, vector<unsigned> &, vector<int> &, vector<lp::numeric_pair<lp::mpq> >&, vector<lp::mpq>&, lp::lp_settings&, const column_namer&, const vector<lp::column_type >&,
+template lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::lp_core_solver_base(lp::static_matrix<lp::mpq, lp::numeric_pair<lp::mpq> >&, 
+    // vector<lp::numeric_pair<lp::mpq> >&, 
+    vector<unsigned int >&, vector<unsigned> &, vector<int> &, vector<lp::numeric_pair<lp::mpq> >&, vector<lp::mpq>&, lp::lp_settings&, const column_namer&, const vector<lp::column_type >&,
                                                                                                    const vector<lp::numeric_pair<lp::mpq> >&,
                                                                                                    const vector<lp::numeric_pair<lp::mpq> >&);
 template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::print_statistics_with_cost_and_check_that_the_time_is_over(lp::numeric_pair<lp::mpq>, std::ostream&);
@@ -91,7 +93,7 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::upda
 template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::add_delta_to_entering(unsigned int, const lp::numeric_pair<lp::mpq>&);
 template lp::lp_core_solver_base<lp::mpq, lp::mpq>::lp_core_solver_base(
                                                                               lp::static_matrix<lp::mpq, lp::mpq>&,
-                                                                              vector<lp::mpq>&,
+                                                                              //vector<lp::mpq>&,
                                                                               vector<unsigned int >&,
                                                                               vector<unsigned> &, vector<int> &,
                                                                               vector<lp::mpq>&,

src/math/lp/lp_core_solver_base.h

@@ -67,7 +67,7 @@ public:
     indexed_vector<T>     m_pivot_row_of_B_1;  // the pivot row of the reverse of B
     indexed_vector<T>     m_pivot_row; // this is the real pivot row of the simplex tableu
     static_matrix<T, X> & m_A; // the matrix A
-    vector<X> &           m_b; // the right side
+    // vector<X> const &           m_b; // the right side
     vector<unsigned> &    m_basis;
     vector<unsigned>&     m_nbasis;
     vector<int>&          m_basis_heading;
@@ -118,7 +118,7 @@ public:
     unsigned m_n() const { return m_A.column_count(); } // the number of columns in the matrix m_A
 
     lp_core_solver_base(static_matrix<T, X> & A,
-                        vector<X> & b, // the right side vector
+                        //vector<X> & b, // the right side vector
                         vector<unsigned> & basis,
                         vector<unsigned> & nbasis,
                         vector<int> & heading,
@@ -213,7 +213,6 @@ public:
         return !below_bound(x, bound) && !above_bound(x, bound);
     }
 
-
     bool need_to_pivot_to_basis_tableau() const {
         unsigned m = m_A.row_count();
         for (unsigned i = 0; i < m; i++) {
@@ -284,7 +283,6 @@ public:
         return below_bound(m_x[p], m_upper_bounds[p]);
     }
 
-
     bool x_above_upper_bound(unsigned p) const {
         return above_bound(m_x[p], m_upper_bounds[p]);
     }
@@ -310,7 +308,6 @@ public:
 
     bool d_is_not_positive(unsigned j) const;
 
-
     bool time_is_over();
 
     void rs_minus_Anx(vector<X> & rs);
@@ -351,8 +348,6 @@ public:
 
     std::string column_name(unsigned column) const;
 
-    void copy_right_side(vector<X> & rs);
-
     void add_delta_to_xB(vector<X> & del);
 
     void find_error_in_BxB(vector<X>& rs);
@@ -367,8 +362,6 @@ public:
         return ret;
     }
     
-    
-    
     bool snap_column_to_bound(unsigned j) {
         switch (m_column_types[j]) {
         case column_type::fixed:

src/math/lp/lp_core_solver_base_def.h

@@ -28,7 +28,7 @@ namespace lp {
 
 template <typename T, typename X> lp_core_solver_base<T, X>::
 lp_core_solver_base(static_matrix<T, X> & A,
-                    vector<X> & b, // the right side vector
+                    // vector<X> & b, // the right side vector
                     vector<unsigned> & basis,
                     vector<unsigned> & nbasis,
                     vector<int> & heading,
@@ -47,7 +47,7 @@ lp_core_solver_base(static_matrix<T, X> & A,
     m_pivot_row_of_B_1(A.row_count()),
     m_pivot_row(A.column_count()),
     m_A(A),
-    m_b(b),
+    // m_b(b),
     m_basis(basis),
     m_nbasis(nbasis),
     m_basis_heading(heading),
@@ -220,7 +220,7 @@ A_mult_x_is_off() const {
     lp_assert(m_x.size() == m_A.column_count());
     if (numeric_traits<T>::precise()) {
 		for (unsigned i = 0; i < m_m(); i++) {
-            X delta = m_b[i] - m_A.dot_product_with_row(i, m_x);
+            X delta = /*m_b[i] */- m_A.dot_product_with_row(i, m_x);
             if (delta != numeric_traits<X>::zero()) {
                 return true;
             }
@@ -230,8 +230,8 @@ A_mult_x_is_off() const {
     T feps = convert_struct<T, double>::convert(m_settings.refactor_tolerance);
     X one = convert_struct<X, double>::convert(1.0);
     for (unsigned i = 0; i < m_m(); i++) {
-        X delta = abs(m_b[i] - m_A.dot_product_with_row(i, m_x));
+        X delta = abs(/*m_b[i] -*/ m_A.dot_product_with_row(i, m_x));
-        X eps = feps * (one + T(0.1) * abs(m_b[i]));
+        auto eps = feps /* * (one + T(0.1) * abs(m_b[i])) */;
 
         if (delta > eps) {
 #if 0
@@ -263,8 +263,8 @@ A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
     T feps = convert_struct<T, double>::convert(m_settings.refactor_tolerance);
     X one = convert_struct<X, double>::convert(1.0);
     for (unsigned i : index) {
-        X delta = abs(m_b[i] - m_A.dot_product_with_row(i, m_x));
+        X delta = abs(/*m_b[i] -*/ m_A.dot_product_with_row(i, m_x));
-        X eps = feps * (one + T(0.1) * abs(m_b[i]));
+        auto eps = feps /* *(one + T(0.1) * abs(m_b[i])) */;
 
         if (delta > eps) {
 #if 0
@@ -400,7 +400,8 @@ column_is_dual_feasible(unsigned j) const {
     case column_type::lower_bound:
         return x_is_at_lower_bound(j) && d_is_not_negative(j);
     case column_type::upper_bound:
-        lp_assert(false); // impossible case
+        UNREACHABLE();
+        break;
     case column_type::free_column:
         return numeric_traits<X>::is_zero(m_d[j]);
     default:
@@ -441,7 +442,7 @@ template <typename T, typename X> void lp_core_solver_base<T, X>::
 rs_minus_Anx(vector<X> & rs) {
     unsigned row = m_m();
     while (row--) {
-        auto &rsv = rs[row] = m_b[row];
+        auto& rsv = rs[row] = zero_of_type<X>(); //m_b[row];
         for (auto & it : m_A.m_rows[row]) {
             unsigned j = it.var();
             if (m_basis_heading[j] < 0) {
@@ -454,8 +455,7 @@ rs_minus_Anx(vector<X> & rs) {
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
 find_x_by_solving() {
     solve_Ax_eq_b();
-    bool ret=  !A_mult_x_is_off();
+    return !A_mult_x_is_off();
-    return ret;
 }
 
 template <typename T, typename X> bool lp_core_solver_base<T, X>::column_is_feasible(unsigned j) const {
@@ -463,28 +463,12 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::column_is_feas
     switch (this->m_column_types[j]) {
     case column_type::fixed:
     case column_type::boxed:
-        if (this->above_bound(x, this->m_upper_bounds[j])) {
+        return !this->above_bound(x, this->m_upper_bounds[j]) && 
-            return false;
+               !this->below_bound(x, this->m_lower_bounds[j]);
-        } else if (this->below_bound(x, this->m_lower_bounds[j])) {
-            return false;
-        } else {
-            return true;
-        }
-        break;
     case column_type::lower_bound:
-        if (this->below_bound(x, this->m_lower_bounds[j])) {
+        return !this->below_bound(x, this->m_lower_bounds[j]);
-            return false;
-        } else {
-            return true;
-        }
-        break;
     case column_type::upper_bound:
-        if (this->above_bound(x, this->m_upper_bounds[j])) {
+        return !this->above_bound(x, this->m_upper_bounds[j]);
-            return false;
-        } else {
-            return true;
-        }
-        break;
     case column_type::free_column:
         return true;
         break;
@@ -598,7 +582,7 @@ divide_row_by_pivot(unsigned pivot_row, unsigned pivot_col) {
     if (is_zero(coeff)) 
         return false;
     
-    this->m_b[pivot_row] /= coeff;
+    // this->m_b[pivot_row] /= coeff;
     for (unsigned j = 0; j < size; j++) {
         auto & c = row[j];
         if (c.var() != pivot_col) {
@@ -662,58 +646,51 @@ basis_has_no_doubles() const {
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
 non_basis_has_no_doubles() const {
     std::set<int> bm;
-    for (auto j : m_nbasis) {
+    for (auto j : m_nbasis) 
         bm.insert(j);    
-    }
     return bm.size() == m_nbasis.size();
 }
 
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
 basis_is_correctly_represented_in_heading() const {
-    for (unsigned i = 0; i < m_m(); i++) {
+    for (unsigned i = 0; i < m_m(); i++) 
         if (m_basis_heading[m_basis[i]] != static_cast<int>(i))
             return false;    
-    }
     return true;
 }
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
 non_basis_is_correctly_represented_in_heading() const {
-    for (unsigned i = 0; i < m_nbasis.size(); i++) {
+    for (unsigned i = 0; i < m_nbasis.size(); i++) 
         if (m_basis_heading[m_nbasis[i]] !=  - static_cast<int>(i) - 1)
             return false;
-    }
+    
-    for (unsigned j = 0; j < m_A.column_count(); j++) {
+    for (unsigned j = 0; j < m_A.column_count(); j++) 
-        if (m_basis_heading[j] >= 0) {
+        if (m_basis_heading[j] >= 0)
             lp_assert(static_cast<unsigned>(m_basis_heading[j]) < m_A.row_count() && m_basis[m_basis_heading[j]] == j);
-        }
+        
-    }
     return true;
 }
 
 template <typename T, typename X> bool lp_core_solver_base<T, X>::
     basis_heading_is_correct() const {
-    if ( m_A.column_count() > 10 ) { // for the performance reason
+    if ( m_A.column_count() > 10 )  // for the performance reason
         return true;
-    }
+    
     lp_assert(m_basis_heading.size() == m_A.column_count());
     lp_assert(m_basis.size() == m_A.row_count());
     lp_assert(m_nbasis.size() <= m_A.column_count() - m_A.row_count()); // for the dual the size of non basis can be smaller
-    if (!basis_has_no_doubles()) {
-        return false;
-    }
 
-    if (!non_basis_has_no_doubles()) {
+    if (!basis_has_no_doubles()) 
         return false;
-    }
     
-    if (!basis_is_correctly_represented_in_heading()) {
+    if (!non_basis_has_no_doubles()) 
         return false;
-    }
     
-    if (!non_basis_is_correctly_represented_in_heading()) {
+    if (!basis_is_correctly_represented_in_heading()) 
         return false;    
-    }
 
+    if (!non_basis_is_correctly_represented_in_heading()) 
+        return false;
     
     return true;
 }
@@ -781,14 +758,6 @@ column_name(unsigned column) const {
     return m_column_names.get_variable_name(column);
 }
 
-template <typename T, typename X> void lp_core_solver_base<T, X>::
-copy_right_side(vector<X> & rs) {
-    unsigned i = m_m();
-    while (i --) {
-        rs[i] = m_b[i];
-    }
-}
-
 template <typename T, typename X> void lp_core_solver_base<T, X>::
 add_delta_to_xB(vector<X> & del) {
     unsigned i = m_m();
@@ -819,7 +788,8 @@ solve_Ax_eq_b() {
         rs_minus_Anx(rs);
         m_factorization->solve_By(rs);
         copy_rs_to_xB(rs);
-    } else {
+    } 
+    else {
         vector<X> rs(m_m());
         rs_minus_Anx(rs);
         vector<X> rrs = rs; // another copy of rs

src/math/lp/lp_dual_core_solver.h

@@ -61,7 +61,7 @@ public:
                         lp_settings & settings,
                         const column_namer & column_names):
         lp_core_solver_base<T, X>(A,
-                                  b,
+                                  // b,
                                   basis,
                                   nbasis,
                                   heading,

src/math/lp/lp_primal_core_solver.h

@@ -948,7 +948,7 @@ public:
                           const vector<X> & upper_bound_values,
                           lp_settings & settings,
                           const column_namer& column_names):
-        lp_core_solver_base<T, X>(A, b,
+        lp_core_solver_base<T, X>(A, // b,
                                   basis,
                                   nbasis,
                                   heading,
@@ -983,7 +983,7 @@ public:
                           const vector<X> & upper_bound_values,
                           lp_settings & settings,
                           const column_namer& column_names):
-        lp_core_solver_base<T, X>(A, b,
+        lp_core_solver_base<T, X>(A, // b,
                                   basis,
                                   nbasis,
                                   heading,