From 25f103db1abc5d76412f85ddeb74d24bcb453f5e Mon Sep 17 00:00:00 2001
From: Lev Nachmanson <levnach@hotmail.com>
Date: Sat, 4 Mar 2023 14:58:49 -0800
Subject: [PATCH] rm_lp
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
---
src/math/lp/core_solver_pretty_printer.h | 1 -
src/math/lp/core_solver_pretty_printer_def.h | 67 +----
src/math/lp/int_solver.cpp | 1 -
src/math/lp/lar_core_solver.h | 68 +----
src/math/lp/lar_core_solver_def.h | 37 +--
src/math/lp/lar_solver.cpp | 57 ++--
src/math/lp/lar_solver.h | 22 +-
src/math/lp/lp_core_solver_base.cpp | 14 +-
src/math/lp/lp_core_solver_base.h | 12 +-
src/math/lp/lp_core_solver_base_def.h | 68 +----
src/math/lp/lp_dual_core_solver_def.h | 22 +-
src/math/lp/lp_primal_core_solver.h | 3 +-
src/math/lp/lp_primal_core_solver_def.h | 244 +-----------------
.../lp/lp_primal_core_solver_tableau_def.h | 34 +--
src/math/lp/lp_settings.h | 6 +-
src/math/lp/lu.cpp | 3 -
src/math/lp/lu.h | 2 -
src/math/lp/lu_def.h | 15 +-
src/test/lp/lp.cpp | 30 +--
19 files changed, 65 insertions(+), 641 deletions(-)
diff --git a/src/math/lp/core_solver_pretty_printer.h b/src/math/lp/core_solver_pretty_printer.h
index 3c0563c32bb..353212dcd1a 100644
--- a/src/math/lp/core_solver_pretty_printer.h
+++ b/src/math/lp/core_solver_pretty_printer.h
@@ -110,7 +110,6 @@ class core_solver_pretty_printer {
return T_to_string(m_exact_column_norms[col]);
}
- void print_exact_norms();
void print_approx_norms();
diff --git a/src/math/lp/core_solver_pretty_printer_def.h b/src/math/lp/core_solver_pretty_printer_def.h
index 23417b691b6..6971061832b 100644
--- a/src/math/lp/core_solver_pretty_printer_def.h
+++ b/src/math/lp/core_solver_pretty_printer_def.h
@@ -59,22 +59,13 @@ core_solver_pretty_printer<T, X>::core_solver_pretty_printer(const lp_core_solve
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_costs() {
- if (!m_core_solver.use_tableau()) {
- vector<T> local_y(m_core_solver.m_m());
- m_core_solver.solve_yB(local_y);
- for (unsigned i = 0; i < ncols(); i++) {
- if (m_core_solver.m_basis_heading[i] < 0) {
- T t = m_core_solver.m_costs[i] - m_core_solver.m_A.dot_product_with_column(local_y, i);
- set_coeff(m_costs, m_cost_signs, i, t, m_core_solver.column_name(i));
- }
- }
- } else {
+
for (unsigned i = 0; i < ncols(); i++) {
if (m_core_solver.m_basis_heading[i] < 0) {
set_coeff(m_costs, m_cost_signs, i, m_core_solver.m_d[i], m_core_solver.column_name(i));
}
}
- }
+
}
template <typename T, typename X> core_solver_pretty_printer<T, X>::~core_solver_pretty_printer() {
@@ -97,7 +88,7 @@ template <typename T, typename X> T core_solver_pretty_printer<T, X>::current_co
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_A_and_signs() {
- if (numeric_traits<T>::precise() && m_core_solver.m_settings.use_tableau()) {
+ if (numeric_traits<T>::precise() ) {
for (unsigned column = 0; column < ncols(); column++) {
vector<T> t(nrows(), zero_of_type<T>());
for (const auto & c : m_core_solver.m_A.m_columns[column]){
@@ -125,23 +116,7 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_m_
m_rs[row] += t[row] * m_core_solver.m_x[column];
}
}
- } else {
- for (unsigned column = 0; column < ncols(); column++) {
- m_core_solver.solve_Bd(column, m_ed_buff, m_w_buff); // puts the result into m_core_solver.m_ed
- string name = m_core_solver.column_name(column);
- for (unsigned row = 0; row < nrows(); row ++) {
- set_coeff(
- m_A[row],
- m_signs[row],
- column,
- m_ed_buff[row],
- name);
- m_rs[row] += m_ed_buff[row] * m_core_solver.m_x[column];
- }
- if (!m_core_solver.use_tableau())
- m_exact_column_norms.push_back(current_column_norm() + T(1)); // a conversion missing 1 -> T
- }
- }
+ }
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::init_column_widths() {
@@ -190,11 +165,7 @@ template <typename T, typename X> unsigned core_solver_pretty_printer<T, X>:: ge
w = cellw;
}
}
- if (!m_core_solver.use_tableau()) {
- w = std::max(w, (unsigned)T_to_string(m_exact_column_norms[column]).size());
- if (!m_core_solver.m_column_norms.empty())
- w = std::max(w, (unsigned)T_to_string(m_core_solver.m_column_norms[column]).size());
- }
+
return w;
}
@@ -315,41 +286,15 @@ template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_u
m_out << std::endl;
}
-template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_exact_norms() {
- if (m_core_solver.use_tableau()) return;
- int blanks = m_title_width + 1 - static_cast<int>(m_exact_norm_title.size());
- m_out << m_exact_norm_title;
- print_blanks_local(blanks, m_out);
- for (unsigned i = 0; i < ncols(); i++) {
- string s = get_exact_column_norm_string(i);
- int blanks = m_column_widths[i] - static_cast<int>(s.size());
- print_blanks_local(blanks, m_out);
- m_out << s << " ";
- }
- m_out << std::endl;
-}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print_approx_norms() {
- if (m_core_solver.use_tableau()) return;
- int blanks = m_title_width + 1 - static_cast<int>(m_approx_norm_title.size());
- m_out << m_approx_norm_title;
- print_blanks_local(blanks, m_out);
- for (unsigned i = 0; i < ncols(); i++) {
- string s = T_to_string(m_core_solver.m_column_norms[i]);
- int blanks = m_column_widths[i] - static_cast<int>(s.size());
- print_blanks_local(blanks, m_out);
- m_out << s << " ";
- }
- m_out << std::endl;
+ return;
}
template <typename T, typename X> void core_solver_pretty_printer<T, X>::print() {
for (unsigned i = 0; i < nrows(); i++) {
print_row(i);
}
- print_exact_norms();
- if (!m_core_solver.m_column_norms.empty())
- print_approx_norms();
m_out << std::endl;
if (m_core_solver.inf_set().size()) {
m_out << "inf columns: ";
diff --git a/src/math/lp/int_solver.cpp b/src/math/lp/int_solver.cpp
index 3aaf8f29eb7..e338e222a27 100644
--- a/src/math/lp/int_solver.cpp
+++ b/src/math/lp/int_solver.cpp
@@ -344,7 +344,6 @@ bool int_solver::get_freedom_interval_for_column(unsigned j, bool & inf_l, impq
set_upper(u, inf_u, upper_bound(j) - xj);
- lp_assert(settings().use_tableau());
const auto & A = lra.A_r();
TRACE("random_update", tout << "m = " << m << "\n";);
diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h
index b6c97513919..9678edd6bed 100644
--- a/src/math/lp/lar_core_solver.h
+++ b/src/math/lp/lar_core_solver.h
@@ -80,8 +80,7 @@ class lar_core_solver {
column_type get_column_type(unsigned j) { return m_column_types[j];}
- void calculate_pivot_row(unsigned i);
-
+
void print_pivot_row(std::ostream & out, unsigned row_index) const {
for (unsigned j : m_r_solver.m_pivot_row.m_index) {
if (numeric_traits<mpq>::is_pos(m_r_solver.m_pivot_row.m_data[j]))
@@ -138,18 +137,11 @@ class lar_core_solver {
void fill_not_improvable_zero_sum();
void pop_basis(unsigned k) {
- if (!settings().use_tableau()) {
- m_r_pushed_basis.pop(k);
- m_r_basis = m_r_pushed_basis();
- m_r_solver.init_basis_heading_and_non_basic_columns_vector();
- m_d_pushed_basis.pop(k);
- m_d_basis = m_d_pushed_basis();
- m_d_solver.init_basis_heading_and_non_basic_columns_vector();
- } else {
+
m_d_basis = m_r_basis;
m_d_nbasis = m_r_nbasis;
m_d_heading = m_r_heading;
- }
+
}
void push() {
@@ -160,19 +152,11 @@ class lar_core_solver {
m_stacked_simplex_strategy.push();
m_column_types.push();
// rational
- if (!settings().use_tableau())
- m_r_A.push();
m_r_lower_bounds.push();
m_r_upper_bounds.push();
- if (!settings().use_tableau()) {
- push_vector(m_r_pushed_basis, m_r_basis);
- push_vector(m_r_columns_nz, m_r_solver.m_columns_nz);
- push_vector(m_r_rows_nz, m_r_solver.m_rows_nz);
- }
m_d_A.push();
- if (!settings().use_tableau())
- push_vector(m_d_pushed_basis, m_d_basis);
+
}
template <typename K>
@@ -202,8 +186,6 @@ class lar_core_solver {
void pop(unsigned k) {
// rationals
- if (!settings().use_tableau())
- m_r_A.pop(k);
m_r_lower_bounds.pop(k);
m_r_upper_bounds.pop(k);
m_column_types.pop(k);
@@ -213,8 +195,7 @@ class lar_core_solver {
m_r_x.resize(m_r_A.column_count());
m_r_solver.m_costs.resize(m_r_A.column_count());
m_r_solver.m_d.resize(m_r_A.column_count());
- if(!settings().use_tableau())
- pop_markowitz_counts(k);
+
m_d_A.pop(k);
// doubles
delete m_d_solver.m_factorization;
@@ -454,7 +435,6 @@ class lar_core_solver {
void solve_on_signature_tableau(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
r_basis_is_OK();
- lp_assert(settings().use_tableau());
bool r = catch_up_in_lu_tableau(changes_of_basis, m_d_solver.m_basis_heading);
if (!r) { // it is the case where m_d_solver gives a degenerated basis
@@ -553,8 +533,7 @@ class lar_core_solver {
bool r_basis_is_OK() const {
#ifdef Z3DEBUG
- if (!m_r_solver.m_settings.use_tableau())
- return true;
+
for (unsigned j : m_r_solver.m_basis) {
lp_assert(m_r_solver.m_A.m_columns[j].size() == 1);
}
@@ -568,40 +547,7 @@ class lar_core_solver {
return true;
}
- void solve_on_signature(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) {
- SASSERT(!settings().use_tableau());
- if (m_r_solver.m_factorization == nullptr) {
- for (unsigned j = 0; j < changes_of_basis.size(); j+=2) {
- unsigned entering = changes_of_basis[j];
- unsigned leaving = changes_of_basis[j + 1];
- m_r_solver.change_basis_unconditionally(entering, leaving);
- }
- init_factorization(m_r_solver.m_factorization, m_r_A, m_r_basis, settings());
- } else {
- catch_up_in_lu(changes_of_basis, m_d_solver.m_basis_heading, m_r_solver);
- }
-
- if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back
- catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver);
- m_r_solver.find_feasible_solution();
- m_d_basis = m_r_basis;
- m_d_heading = m_r_heading;
- m_d_nbasis = m_r_nbasis;
- delete m_d_solver.m_factorization;
- m_d_solver.m_factorization = nullptr;
- } else {
- prepare_solver_x_with_signature(signature, m_r_solver);
- m_r_solver.start_tracing_basis_changes();
- m_r_solver.find_feasible_solution();
- if (settings().get_cancel_flag())
- return;
- m_r_solver.stop_tracing_basis_changes();
- // and now catch up in the double solver
- lp_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2);
- catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver);
- }
- }
-
+
void create_double_matrix(static_matrix<double, double> & A) {
for (unsigned i = 0; i < m_r_A.row_count(); i++) {
auto & row = m_r_A.m_rows[i];
diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h
index 75fff64fd71..182029e3e98 100644
--- a/src/math/lp/lar_core_solver_def.h
+++ b/src/math/lp/lar_core_solver_def.h
@@ -46,23 +46,9 @@ 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());
- m_r_solver.m_ed.resize(m_r_solver.m_m());
- m_r_solver.m_pivot_row.resize(m_r_solver.m_n());
- m_r_solver.m_pivot_row_of_B_1.resize(m_r_solver.m_m());
- m_r_solver.m_w.resize(m_r_solver.m_m());
- m_r_solver.m_y.resize(m_r_solver.m_m());
- m_r_solver.m_rows_nz.resize(m_r_solver.m_m(), 0);
- m_r_solver.m_columns_nz.resize(m_r_solver.m_n(), 0);
- init_column_row_nz_for_r_solver();
- }
-
+
// 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)
@@ -142,7 +128,7 @@ void lar_core_solver::solve() {
}
++settings().stats().m_need_to_solve_inf;
CASSERT("A_off", !m_r_solver.A_mult_x_is_off());
- lp_assert((!settings().use_tableau()) || r_basis_is_OK());
+ lp_assert( r_basis_is_OK());
if (need_to_presolve_with_double_solver()) {
TRACE("lar_solver", tout << "presolving\n";);
prefix_d();
@@ -152,26 +138,17 @@ void lar_core_solver::solve() {
m_r_solver.set_status(lp_status::TIME_EXHAUSTED);
return;
}
- if (settings().use_tableau())
- solve_on_signature_tableau(solution_signature, changes_of_basis);
- else
- solve_on_signature(solution_signature, changes_of_basis);
-
- lp_assert(!settings().use_tableau() || r_basis_is_OK());
+ solve_on_signature_tableau(solution_signature, changes_of_basis);
+
+ lp_assert( r_basis_is_OK());
} else {
- if (!settings().use_tableau()) {
- TRACE("lar_solver", tout << "no tablau\n";);
- bool snapped = m_r_solver.snap_non_basic_x_to_bound();
- lp_assert(m_r_solver.non_basic_columns_are_set_correctly());
- if (snapped)
- m_r_solver.solve_Ax_eq_b();
- }
+
if (m_r_solver.m_look_for_feasible_solution_only) //todo : should it be set?
m_r_solver.find_feasible_solution();
else {
m_r_solver.solve();
}
- lp_assert(!settings().use_tableau() || r_basis_is_OK());
+ lp_assert(r_basis_is_OK());
}
switch (m_r_solver.get_status())
{
diff --git a/src/math/lp/lar_solver.cpp b/src/math/lp/lar_solver.cpp
index 3339b55a5b7..2e9541525d4 100644
--- a/src/math/lp/lar_solver.cpp
+++ b/src/math/lp/lar_solver.cpp
@@ -261,8 +261,7 @@ namespace lp {
m_crossed_bounds_column.pop(k);
unsigned n = m_columns_to_ul_pairs.peek_size(k);
m_var_register.shrink(n);
- if (m_settings.use_tableau())
- pop_tableau();
+ pop_tableau();
lp_assert(A_r().column_count() == n);
TRACE("lar_solver_details",
for (unsigned j = 0; j < n; j++) {
@@ -284,7 +283,7 @@ namespace lp {
clean_popped_elements(m, m_rows_with_changed_bounds);
clean_inf_set_of_r_solver_after_pop();
lp_assert(m_settings.simplex_strategy() == simplex_strategy_enum::undecided ||
- (!use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
+ ( m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau()));
m_constraints.pop(k);
@@ -299,7 +298,7 @@ namespace lp {
m_simplex_strategy.pop(k);
m_settings.set_simplex_strategy(m_simplex_strategy);
lp_assert(sizes_are_correct());
- lp_assert((!m_settings.use_tableau()) || m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
+ lp_assert(m_mpq_lar_core_solver.m_r_solver.reduced_costs_are_correct_tableau());
m_usage_in_terms.pop(k);
set_status(lp_status::UNKNOWN);
}
@@ -630,15 +629,7 @@ namespace lp {
}
void lar_solver::detect_rows_of_bound_change_column_for_nbasic_column(unsigned j) {
- if (A_r().row_count() != m_column_buffer.data_size())
- m_column_buffer.resize(A_r().row_count());
- else
- m_column_buffer.clear();
- lp_assert(m_column_buffer.size() == 0 && m_column_buffer.is_OK());
-
- m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer);
- for (unsigned i : m_column_buffer.m_index)
- insert_row_with_changed_bounds(i);
+ lp_assert(false);
}
@@ -648,8 +639,7 @@ namespace lp {
insert_row_with_changed_bounds(rc.var());
}
- bool lar_solver::use_tableau() const { return m_settings.use_tableau(); }
-
+
bool lar_solver::use_tableau_costs() const {
return m_settings.simplex_strategy() == simplex_strategy_enum::tableau_costs;
}
@@ -692,7 +682,7 @@ namespace lp {
}
void lar_solver::change_basic_columns_dependend_on_a_given_nb_column(unsigned j, const numeric_pair<mpq>& delta) {
- if (use_tableau()) {
+ {
for (const auto& c : A_r().m_columns[j]) {
unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
if (tableau_with_costs()) {
@@ -705,15 +695,7 @@ namespace lp {
}
}
- else {
- m_column_buffer.clear();
- m_column_buffer.resize(A_r().row_count());
- m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer);
- for (unsigned i : m_column_buffer.m_index) {
- unsigned bj = m_mpq_lar_core_solver.m_r_basis[i];
- m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, -m_column_buffer[i] * delta);
- }
- }
+
}
void lar_solver::update_x_and_inf_costs_for_column_with_changed_bounds(unsigned j) {
@@ -742,10 +724,8 @@ namespace lp {
return;
}
- if (use_tableau())
- detect_rows_of_bound_change_column_for_nbasic_column_tableau(j);
- else
- detect_rows_of_bound_change_column_for_nbasic_column(j);
+ detect_rows_of_bound_change_column_for_nbasic_column_tableau(j);
+
}
void lar_solver::detect_rows_with_changed_bounds() {
@@ -773,8 +753,6 @@ namespace lp {
void lar_solver::solve_with_core_solver() {
- if (!use_tableau())
- add_last_rows_to_lu(m_mpq_lar_core_solver.m_r_solver);
if (m_mpq_lar_core_solver.need_to_presolve_with_double_solver()) {
add_last_rows_to_lu(m_mpq_lar_core_solver.m_d_solver);
}
@@ -782,10 +760,7 @@ namespace lp {
if (costs_are_used()) {
m_basic_columns_with_changed_cost.resize(m_mpq_lar_core_solver.m_r_x.size());
}
- if (use_tableau())
- update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();
- else
- update_x_and_inf_costs_for_columns_with_changed_bounds();
+ update_x_and_inf_costs_for_columns_with_changed_bounds_tableau();
m_mpq_lar_core_solver.solve();
set_status(m_mpq_lar_core_solver.m_r_solver.get_status());
lp_assert(((stats().m_make_feasible% 100) != 0) || m_status != lp_status::OPTIMAL || all_constraints_hold());
@@ -1753,7 +1728,7 @@ namespace lp {
SASSERT(m_terms.size() == m_term_register.size());
unsigned adjusted_term_index = m_terms.size() - 1;
var_index ret = tv::mask_term(adjusted_term_index);
- if (use_tableau() && !coeffs.empty()) {
+ if ( !coeffs.empty()) {
add_row_from_term_no_constraint(m_terms.back(), ret);
if (m_settings.bound_propagation())
insert_row_with_changed_bounds(A_r().row_count() - 1);
@@ -1774,14 +1749,12 @@ namespace lp {
ul_pair ul(true); // to mark this column as associated_with_row
m_columns_to_ul_pairs.push_back(ul);
add_basic_var_to_core_fields();
- if (use_tableau()) {
- A_r().fill_last_row_with_pivoting(*term,
+
+ A_r().fill_last_row_with_pivoting(*term,
j,
m_mpq_lar_core_solver.m_r_solver.m_basis_heading);
- }
- else {
- fill_last_row_of_A_r(A_r(), term);
- }
+
+
m_mpq_lar_core_solver.m_r_solver.update_x(j, get_basic_var_value_from_row(A_r().row_count() - 1));
for (lar_term::ival c : *term) {
unsigned j = c.column();
diff --git a/src/math/lp/lar_solver.h b/src/math/lp/lar_solver.h
index c1deb2fb1c8..51365936819 100644
--- a/src/math/lp/lar_solver.h
+++ b/src/math/lp/lar_solver.h
@@ -177,7 +177,6 @@ class lar_solver : public column_namer {
if (A_r().m_rows[row_index].size() > settings().max_row_length_for_bound_propagation
|| row_has_a_big_num(row_index))
return;
- lp_assert(use_tableau());
bound_analyzer_on_row<row_strip<mpq>, lp_bound_propagator<T>>::analyze_row(A_r().m_rows[row_index],
null_ci,
@@ -190,7 +189,6 @@ class lar_solver : public column_namer {
void substitute_basis_var_in_terms_for_row(unsigned i);
template <typename T>
void calculate_implied_bounds_for_row(unsigned i, lp_bound_propagator<T> & bp) {
- SASSERT(use_tableau());
analyze_new_bounds_on_row_tableau(i, bp);
}
static void clean_popped_elements(unsigned n, u_set& set);
@@ -210,7 +208,6 @@ class lar_solver : public column_namer {
vector<std::pair<mpq, var_index>> &left_side) const;
void detect_rows_of_bound_change_column_for_nbasic_column(unsigned j);
void detect_rows_of_bound_change_column_for_nbasic_column_tableau(unsigned j);
- bool use_tableau() const;
bool use_tableau_costs() const;
void detect_rows_of_column_with_bound_change(unsigned j);
void adjust_x_of_column(unsigned j);
@@ -376,7 +373,6 @@ class lar_solver : public column_namer {
void mark_rows_for_bound_prop(lpvar j);
template <typename T>
void propagate_bounds_for_touched_rows(lp_bound_propagator<T> & bp) {
- SASSERT(use_tableau());
for (unsigned i : m_rows_with_changed_bounds) {
calculate_implied_bounds_for_row(i, bp);
if (settings().get_cancel_flag())
@@ -429,8 +425,8 @@ class lar_solver : public column_namer {
void change_basic_columns_dependend_on_a_given_nb_column_report(unsigned j,
const numeric_pair<mpq> & delta,
const ChangeReport& after) {
- if (use_tableau()) {
- for (const auto & c : A_r().m_columns[j]) {
+
+ for (const auto & c : A_r().m_columns[j]) {
unsigned bj = m_mpq_lar_core_solver.m_r_basis[c.var()];
if (tableau_with_costs()) {
m_basic_columns_with_changed_cost.insert(bj);
@@ -440,20 +436,8 @@ class lar_solver : public column_namer {
TRACE("change_x_del",
tout << "changed basis column " << bj << ", it is " <<
( m_mpq_lar_core_solver.m_r_solver.column_is_feasible(bj)? "feas":"inf") << std::endl;);
-
-
- }
- } else {
- NOT_IMPLEMENTED_YET();
- m_column_buffer.clear();
- m_column_buffer.resize(A_r().row_count());
- m_mpq_lar_core_solver.m_r_solver.solve_Bd(j, m_column_buffer);
- for (unsigned i : m_column_buffer.m_index) {
- unsigned bj = m_mpq_lar_core_solver.m_r_basis[i];
- m_mpq_lar_core_solver.m_r_solver.add_delta_to_x_and_track_feasibility(bj, -m_column_buffer[i] * delta);
}
- }
- }
+ }
template <typename ChangeReport>
void set_value_for_nbasic_column_report(unsigned j,
diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp
index c8b1692d272..5dc8fb9e2c2 100644
--- a/src/math/lp/lp_core_solver_base.cpp
+++ b/src/math/lp/lp_core_solver_base.cpp
@@ -26,7 +26,6 @@ Revision History:
template bool lp::lp_core_solver_base<double, double>::A_mult_x_is_off() const;
template bool lp::lp_core_solver_base<double, double>::A_mult_x_is_off_on_index(const vector<unsigned> &) const;
template bool lp::lp_core_solver_base<double, double>::basis_heading_is_correct() const;
-template void lp::lp_core_solver_base<double, double>::calculate_pivot_row_of_B_1(unsigned int);
template void lp::lp_core_solver_base<double, double>::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned);
template bool lp::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int) const;
template void lp::lp_core_solver_base<double, double>::fill_reduced_costs_from_m_y_by_rows();
@@ -52,16 +51,12 @@ template void lp::lp_core_solver_base<double, double>::set_non_basic_x_to_correc
template void lp::lp_core_solver_base<double, double>::snap_xN_to_bounds_and_free_columns_to_zeroes();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::snap_xN_to_bounds_and_free_columns_to_zeroes();
template void lp::lp_core_solver_base<double, double>::solve_Ax_eq_b();
-template void lp::lp_core_solver_base<double, double>::solve_Bd(unsigned int);
-template void lp::lp_core_solver_base<double, double>::solve_Bd(unsigned int, lp::indexed_vector<double>&, lp::indexed_vector<double>&) const;
-template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq>>::solve_Bd(unsigned int, indexed_vector<lp::mpq>&);
template void lp::lp_core_solver_base<double, double>::solve_yB(vector<double >&) const;
template bool lp::lp_core_solver_base<double, double>::update_basis_and_x(int, int, double const&);
template void lp::lp_core_solver_base<double, double>::add_delta_to_entering(unsigned int, const double&);
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::A_mult_x_is_off() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::A_mult_x_is_off_on_index(const vector<unsigned> &) const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::basis_heading_is_correct() const ;
-template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::calculate_pivot_row_of_B_1(unsigned int);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned);
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::column_is_dual_feasible(unsigned int) const;
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::fill_reduced_costs_from_m_y_by_rows();
@@ -71,11 +66,9 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::print_statistics_with_i
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::restore_x(unsigned int, lp::mpq const&);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::set_non_basic_x_to_correct_bounds();
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::solve_Ax_eq_b();
-template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::solve_Bd(unsigned int);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::solve_yB(vector<lp::mpq>&) const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::update_basis_and_x(int, int, lp::mpq const&);
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::add_delta_to_entering(unsigned int, const lp::mpq&);
-template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calculate_pivot_row_of_B_1(unsigned int);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned);
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();
@@ -88,7 +81,7 @@ template lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::lp_core_s
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&);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::snap_xN_to_bounds_and_fill_xB();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_Ax_eq_b();
-template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_Bd(unsigned int);
+
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::update_basis_and_x(int, int, lp::numeric_pair<lp::mpq> const&);
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(
@@ -145,8 +138,7 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::inf_set_is_correct() co
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::infeasibility_costs_are_correct() const;
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq >::infeasibility_costs_are_correct() const;
template bool lp::lp_core_solver_base<double, double >::infeasibility_costs_are_correct() const;
-template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::calculate_pivot_row(unsigned int);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_from_basis(unsigned int);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::remove_from_basis(unsigned int, lp::numeric_pair<lp::mpq> const&);
-template void lp::lp_core_solver_base<rational, rational>::solve_Bd(unsigned int, lp::indexed_vector<rational>&, lp::indexed_vector<rational>&) const;
-template void lp::lp_core_solver_base<rational, lp::numeric_pair<rational> >::solve_Bd(unsigned int, lp::indexed_vector<rational>&, lp::indexed_vector<rational>&) const;
+
+
diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h
index b7010aa54a8..631f687816c 100644
--- a/src/math/lp/lp_core_solver_base.h
+++ b/src/math/lp/lp_core_solver_base.h
@@ -156,11 +156,6 @@ class lp_core_solver_base {
void solve_yB(vector<T> & y) const;
- void solve_Bd(unsigned entering, indexed_vector<T> & d_buff, indexed_vector<T>& w_buff) const;
-
- void solve_Bd(unsigned entering);
-
- void solve_Bd(unsigned entering, indexed_vector<T> & column);
void pretty_print(std::ostream & out);
@@ -184,9 +179,7 @@ class lp_core_solver_base {
bool A_mult_x_is_off() const;
bool A_mult_x_is_off_on_index(const vector<unsigned> & index) const;
- // from page 182 of Istvan Maros's book
- void calculate_pivot_row_of_B_1(unsigned pivot_row);
-
+
void calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned pivot_row);
void add_delta_to_entering(unsigned entering, const X & delta);
@@ -670,7 +663,6 @@ class lp_core_solver_base {
m_settings.simplex_strategy();
}
- bool use_tableau() const { return m_settings.use_tableau(); }
template <typename K>
static void swap(vector<K> &v, unsigned i, unsigned j) {
@@ -760,7 +752,7 @@ class lp_core_solver_base {
return m_iters_with_no_cost_growing;
}
- void calculate_pivot_row(unsigned i);
+
unsigned get_base_column_in_row(unsigned row_index) const {
return m_basis[row_index];
}
diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h
index 93d5f430239..0e834243050 100644
--- a/src/math/lp/lp_core_solver_base_def.h
+++ b/src/math/lp/lp_core_solver_base_def.h
@@ -131,35 +131,9 @@ solve_yB(vector<T> & y) const {
// m_index_of_ed.push_back(i);
// }
// }
-template <typename T, typename X> void lp_core_solver_base<T, X>::solve_Bd(unsigned entering, indexed_vector<T> & column) {
- lp_assert(!m_settings.use_tableau());
- if (m_factorization == nullptr) {
- init_factorization(m_factorization, m_A, m_basis, m_settings);
- }
- m_factorization->solve_Bd_faster(entering, column);
-}
-template <typename T, typename X> void lp_core_solver_base<T, X>::solve_Bd(unsigned , indexed_vector<T>& , indexed_vector<T> &) const {
- NOT_IMPLEMENTED_YET();
-}
-template <typename T, typename X> void lp_core_solver_base<T, X>::
-solve_Bd(unsigned entering) {
- lp_assert(m_ed.is_OK());
- m_factorization->solve_Bd(entering, m_ed, m_w);
- if (this->precise())
- m_columns_nz[entering] = m_ed.m_index.size();
- lp_assert(m_ed.is_OK());
- lp_assert(m_w.is_OK());
-#ifdef Z3DEBUG
- // auto B = get_B(*m_factorization, m_basis);
- // vector<T> a(m_m());
- // m_A.copy_column_to_vector(entering, a);
- // vector<T> cd(m_ed.m_data);
- // B.apply_from_left(cd, m_settings);
- // lp_assert(vectors_are_equal(cd , a));
-#endif
-}
+
template <typename T, typename X> void lp_core_solver_base<T, X>::
pretty_print(std::ostream & out) {
@@ -279,17 +253,6 @@ A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
return false;
}
-// from page 182 of Istvan Maros's book
-template <typename T, typename X> void lp_core_solver_base<T, X>::
-calculate_pivot_row_of_B_1(unsigned pivot_row) {
- lp_assert(! use_tableau());
- lp_assert(m_pivot_row_of_B_1.is_OK());
- m_pivot_row_of_B_1.clear();
- m_pivot_row_of_B_1.set_value(numeric_traits<T>::one(), pivot_row);
- lp_assert(m_pivot_row_of_B_1.is_OK());
- m_factorization->solve_yB_with_error_check_indexed(m_pivot_row_of_B_1, m_basis_heading, m_basis, m_settings);
- lp_assert(m_pivot_row_of_B_1.is_OK());
-}
template <typename T, typename X> void lp_core_solver_base<T, X>::
@@ -316,13 +279,7 @@ calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned pivot_row) {
template <typename T, typename X> void lp_core_solver_base<T, X>::
add_delta_to_entering(unsigned entering, const X& delta) {
m_x[entering] += delta;
- if (!use_tableau())
- for (unsigned i : m_ed.m_index) {
- if (!numeric_traits<X>::precise())
- m_copy_of_xB[i] = m_x[m_basis[i]];
- m_x[m_basis[i]] -= delta * m_ed[i];
- }
- else
+
for (const auto & c : m_A.m_columns[entering]) {
unsigned i = c.var();
m_x[m_basis[i]] -= delta * m_A.get_val(c);
@@ -1000,26 +957,5 @@ lp_core_solver_base<T, X>::infeasibility_cost_is_correct_for_column(unsigned j)
}
}
-template <typename T, typename X>
-void lp_core_solver_base<T, X>::calculate_pivot_row(unsigned i) {
- lp_assert(!use_tableau());
- lp_assert(m_pivot_row.is_OK());
- m_pivot_row_of_B_1.clear();
- m_pivot_row_of_B_1.resize(m_m());
- m_pivot_row.clear();
- m_pivot_row.resize(m_n());
- if (m_settings.use_tableau()) {
- unsigned basic_j = m_basis[i];
- for (auto & c : m_A.m_rows[i]) {
- if (c.var() != basic_j)
- m_pivot_row.set_value(c.coeff(), c.var());
- }
- return;
- }
-
- calculate_pivot_row_of_B_1(i);
- calculate_pivot_row_when_pivot_row_of_B1_is_ready(i);
-}
-
}
diff --git a/src/math/lp/lp_dual_core_solver_def.h b/src/math/lp/lp_dual_core_solver_def.h
index b42d644af34..df70e64f1bc 100644
--- a/src/math/lp/lp_dual_core_solver_def.h
+++ b/src/math/lp/lp_dual_core_solver_def.h
@@ -210,26 +210,8 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::DSE_FTran() {
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_known_p() {
- if (done()) {
- return true;
- }
- this->calculate_pivot_row_of_B_1(m_r);
- this->calculate_pivot_row_when_pivot_row_of_B1_is_ready(m_r);
- if (!ratio_test()) {
- return true;
- }
- calculate_beta_r_precisely();
- 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
- lp_unreachable(); // not implemented yet
- } else {
- lp_assert(pivot_compare_result == 1);
- this->init_lu();
- }
- DSE_FTran();
- return basis_change_and_update();
+
+ return false;
}
template <typename T, typename X> int lp_dual_core_solver<T, X>::define_sign_of_alpha_r() {
diff --git a/src/math/lp/lp_primal_core_solver.h b/src/math/lp/lp_primal_core_solver.h
index a60395ab014..dc6cb2900e8 100644
--- a/src/math/lp/lp_primal_core_solver.h
+++ b/src/math/lp/lp_primal_core_solver.h
@@ -468,7 +468,6 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
}
void update_basis_and_x_tableau_rows(int entering, int leaving, X const & tt) {
- lp_assert(this->use_tableau());
lp_assert(entering != leaving);
update_x_tableau_rows(entering, leaving, tt);
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
@@ -804,7 +803,7 @@ class lp_primal_core_solver:public lp_core_solver_base<T, X> {
return (a > zero_of_type<L>() && m_sign_of_entering_delta > 0) || (a < zero_of_type<L>() && m_sign_of_entering_delta < 0);
}
- void init_reduced_costs();
+
bool lower_bounds_are_set() const override { return true; }
diff --git a/src/math/lp/lp_primal_core_solver_def.h b/src/math/lp/lp_primal_core_solver_def.h
index 3818b589a5a..0a58b0fdf69 100644
--- a/src/math/lp/lp_primal_core_solver_def.h
+++ b/src/math/lp/lp_primal_core_solver_def.h
@@ -33,21 +33,14 @@ namespace lp {
template <typename T, typename X>
void lp_primal_core_solver<T, X>::sort_non_basis_rational() {
lp_assert(numeric_traits<T>::precise());
- if (this->m_settings.use_tableau()) {
+
std::sort(this->m_nbasis.begin(), this->m_nbasis.end(), [this](unsigned a, unsigned b) {
unsigned ca = this->m_A.number_of_non_zeroes_in_column(a);
unsigned cb = this->m_A.number_of_non_zeroes_in_column(b);
if (ca == 0 && cb != 0) return false;
return ca < cb;
});
- } else {
- std::sort(this->m_nbasis.begin(), this->m_nbasis.end(), [this](unsigned a, unsigned b) {
- unsigned ca = this->m_columns_nz[a];
- unsigned cb = this->m_columns_nz[b];
- if (ca == 0 && cb != 0) return false;
- return ca < cb;
- });}
-
+
m_non_basis_list.clear();
// reinit m_basis_heading
for (unsigned j = 0; j < this->m_nbasis.size(); j++) {
@@ -644,25 +637,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::backup_an
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run() {
- this->m_basis_sort_counter = 0; // to initiate the sort of the basis
- // this->set_total_iterations(0);
- this->iters_with_no_cost_growing() = 0;
- init_inf_set();
- if (this->current_x_is_feasible() && this->m_look_for_feasible_solution_only)
- return;
- this->set_using_infeas_costs(false);
- if (this->m_settings.backup_costs)
- backup_and_normalize_costs();
- m_epsilon_of_reduced_cost = numeric_traits<X>::precise()? zero_of_type<T>(): T(1)/T(10000000);
- m_breakpoint_indices_queue.resize(this->m_n());
- init_reduced_costs();
- if (!numeric_traits<X>::precise()) {
- this->m_column_norm_update_counter = 0;
- init_column_norms();
- } else {
- if (this->m_columns_nz.size() != this->m_n())
- init_column_row_non_zeroes();
- }
+
}
@@ -676,166 +651,20 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::calc_work
template <typename T, typename X>
void lp_primal_core_solver<T, X>::advance_on_entering_equal_leaving(int entering, X & t) {
- CASSERT("A_off", !this->A_mult_x_is_off() );
- this->add_delta_to_entering(entering, t * m_sign_of_entering_delta);
- if (this->A_mult_x_is_off_on_index(this->m_ed.m_index) && !this->find_x_by_solving()) {
- this->init_lu();
- if (!this->find_x_by_solving()) {
- this->restore_x(entering, t * m_sign_of_entering_delta);
- this->iters_with_no_cost_growing()++;
- LP_OUT(this->m_settings, "failing in advance_on_entering_equal_leaving for entering = " << entering << std::endl);
- return;
- }
- }
- if (this->using_infeas_costs()) {
- lp_assert(is_zero(this->m_costs[entering]));
- init_infeasibility_costs_for_changed_basis_only();
- }
- if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
- return;
- if (need_to_switch_costs() ||!this->current_x_is_feasible()) {
- init_reduced_costs();
- }
- this->iters_with_no_cost_growing() = 0;
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_entering_and_leaving(int entering, int leaving, X & t) {
- lp_assert(entering >= 0 && m_non_basis_list.back() == static_cast<unsigned>(entering));
- lp_assert(this->using_infeas_costs() || t >= zero_of_type<X>());
- lp_assert(leaving >= 0 && entering >= 0);
- lp_assert(entering != leaving || !is_zero(t)); // otherwise nothing changes
- if (entering == leaving) {
- advance_on_entering_equal_leaving(entering, t);
- return;
- }
- unsigned pivot_row = this->m_basis_heading[leaving];
- this->calculate_pivot_row_of_B_1(pivot_row);
- this->calculate_pivot_row_when_pivot_row_of_B1_is_ready(pivot_row);
-
- 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
- this->set_status(lp_status::UNSTABLE);
- this->iters_with_no_cost_growing()++;
- return;
- } else {
- lp_assert(pivot_compare_result == 1);
- this->init_lu();
- if (this->m_factorization == nullptr || this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::UNSTABLE);
- this->iters_with_no_cost_growing()++;
- return;
- }
- }
- if (!numeric_traits<T>::precise())
- calc_working_vector_beta_for_column_norms();
- if (this->current_x_is_feasible() || !this->m_settings.use_breakpoints_in_feasibility_search) {
- if (m_sign_of_entering_delta == -1)
- t = -t;
- }
- if (!this->update_basis_and_x(entering, leaving, t)) {
- if (this->get_status() == lp_status::FLOATING_POINT_ERROR)
- return;
- if (this->m_look_for_feasible_solution_only) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- return;
- }
- init_reduced_costs();
- return;
- }
-
- if (!is_zero(t)) {
- this->iters_with_no_cost_growing() = 0;
- init_infeasibility_after_update_x_if_inf(leaving);
- }
-
- if (this->current_x_is_feasible()) {
- this->set_status(lp_status::FEASIBLE);
- if (this->m_look_for_feasible_solution_only)
- return;
- }
- if (numeric_traits<X>::precise() == false)
- update_or_init_column_norms(entering, leaving);
-
-
- if (need_to_switch_costs()) {
- init_reduced_costs();
- } else {
- update_reduced_costs_from_pivot_row(entering, leaving);
- }
- lp_assert(!need_to_switch_costs());
- std::list<unsigned>::iterator it = m_non_basis_list.end();
- it--;
- * it = static_cast<unsigned>(leaving);
+
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering_precise(int entering) {
- lp_assert(numeric_traits<T>::precise());
- lp_assert(entering > -1);
- this->solve_Bd(entering);
- X t;
- int leaving = find_leaving_and_t_precise(entering, t);
- if (leaving == -1) {
- TRACE("lar_solver", tout << "non-leaving\n";);
- this->set_status(lp_status::UNBOUNDED);
- return;
- }
- advance_on_entering_and_leaving(entering, leaving, t);
+ lp_assert(false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_on_entering(int entering) {
- if (numeric_traits<T>::precise()) {
- advance_on_entering_precise(entering);
- return;
- }
- lp_assert(entering > -1);
- this->solve_Bd(entering);
- int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
- if (refresh_result) {
- if (this->m_look_for_feasible_solution_only) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- return;
- }
-
- this->init_lu();
- init_reduced_costs();
- if (refresh_result == 2) {
- this->iters_with_no_cost_growing()++;
- return;
- }
- }
- X t;
- int leaving = find_leaving_and_t(entering, t);
- if (leaving == -1){
- if (!this->current_x_is_feasible()) {
- lp_assert(!numeric_traits<T>::precise()); // we cannot have unbounded with inf costs
-
- // if (m_look_for_feasible_solution_only) {
- // this->m_status = INFEASIBLE;
- // return;
- // }
-
-
- if (this->get_status() == lp_status::UNSTABLE) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- return;
- }
- init_infeasibility_costs();
- this->set_status(lp_status::UNSTABLE);
-
- return;
- }
- if (this->get_status() == lp_status::TENTATIVE_UNBOUNDED) {
- this->set_status(lp_status::UNBOUNDED);
- } else {
- this->set_status(lp_status::TENTATIVE_UNBOUNDED);
- }
- TRACE("lar_solver", tout << this->get_status() << "\n";);
- return;
- }
- advance_on_entering_and_leaving(entering, leaving, t);
+ lp_assert(false);
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::push_forward_offset_in_non_basis(unsigned & offset_in_nb) {
@@ -867,7 +696,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_column
// returns the number of iterations
template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve() {
TRACE("lar_solver", tout << "solve " << this->get_status() << "\n";);
- if (numeric_traits<T>::precise() && this->m_settings.use_tableau())
+ if (numeric_traits<T>::precise())
return solve_with_tableau();
init_run();
@@ -893,56 +722,19 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
case lp_status::INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
- if (!numeric_traits<T>::precise()) {
- if(this->m_look_for_feasible_solution_only)
- break;
- this->init_lu();
+ { // precise case
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status (lp_status::FLOATING_POINT_ERROR);
- break;
- }
- init_reduced_costs();
- if (choose_entering_column(1) == -1) {
- decide_on_status_when_cannot_find_entering();
- break;
- }
- this->set_status(lp_status::UNKNOWN);
- } else { // precise case
- if (this->m_look_for_feasible_solution_only) { // todo: keep the reduced costs correct all the time!
- init_reduced_costs();
- if (choose_entering_column(1) == -1) {
- decide_on_status_when_cannot_find_entering();
- break;
- }
- this->set_status(lp_status::UNKNOWN);
- }
}
break;
case lp_status::TENTATIVE_UNBOUNDED:
- this->init_lu();
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- break;
- }
-
- init_reduced_costs();
+ lp_assert(false);
break;
case lp_status::UNBOUNDED:
- if (this->current_x_is_infeasible()) {
- init_reduced_costs();
- this->set_status(lp_status::UNKNOWN);
- }
+ lp_assert(false);
break;
case lp_status::UNSTABLE:
- lp_assert(! (numeric_traits<T>::precise()));
- this->init_lu();
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- break;
- }
- init_reduced_costs();
+ lp_assert(false);
break;
default:
@@ -1292,20 +1084,6 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::print_breakp
print_bound_info_and_x(b->m_j, out);
}
-template <typename T, typename X>
-void lp_primal_core_solver<T, X>::init_reduced_costs() {
- lp_assert(!this->use_tableau());
- if (this->current_x_is_infeasible() && !this->using_infeas_costs()) {
- init_infeasibility_costs();
- } else if (this->current_x_is_feasible() && this->using_infeas_costs()) {
- if (this->m_look_for_feasible_solution_only)
- return;
- this->m_costs = m_costs_backup;
- this->set_using_infeas_costs(false);
- }
-
- this->init_reduced_costs_for_one_iteration();
-}
template <typename T, typename X> void lp_primal_core_solver<T, X>::change_slope_on_breakpoint(unsigned entering, breakpoint<X> * b, T & slope_at_entering) {
if (b->m_j == entering) {
diff --git a/src/math/lp/lp_primal_core_solver_tableau_def.h b/src/math/lp/lp_primal_core_solver_tableau_def.h
index 46297a63e14..fa25694adc7 100644
--- a/src/math/lp/lp_primal_core_solver_tableau_def.h
+++ b/src/math/lp/lp_primal_core_solver_tableau_def.h
@@ -126,22 +126,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
case lp_status::INFEASIBLE:
if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible())
break;
- if (!numeric_traits<T>::precise()) {
- if(this->m_look_for_feasible_solution_only)
- break;
- this->init_lu();
-
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- break;
- }
- init_reduced_costs();
- if (choose_entering_column(1) == -1) {
- decide_on_status_when_cannot_find_entering();
- break;
- }
- this->set_status(lp_status::UNKNOWN);
- } else { // precise case
+ { // precise case
if ((!this->infeasibility_costs_are_correct())) {
init_reduced_costs_tableau(); // forcing recalc
if (choose_entering_column_tableau() == -1) {
@@ -153,13 +138,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
}
break;
case lp_status::TENTATIVE_UNBOUNDED:
- this->init_lu();
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- break;
- }
-
- init_reduced_costs();
+ lp_assert(false);
break;
case lp_status::UNBOUNDED:
if (this->current_x_is_infeasible()) {
@@ -169,13 +148,7 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
break;
case lp_status::UNSTABLE:
- lp_assert(! (numeric_traits<T>::precise()));
- this->init_lu();
- if (this->m_factorization->get_status() != LU_status::OK) {
- this->set_status(lp_status::FLOATING_POINT_ERROR);
- break;
- }
- init_reduced_costs();
+ lp_assert(false);
break;
default:
@@ -348,7 +321,6 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run_tab
template <typename T, typename X> bool lp_primal_core_solver<T, X>::
update_basis_and_x_tableau(int entering, int leaving, X const & tt) {
- lp_assert(this->use_tableau());
lp_assert(entering != leaving);
update_x_tableau(entering, tt);
this->pivot_column_tableau(entering, this->m_basis_heading[leaving]);
diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h
index d6a78564d8f..360ef99bf1d 100644
--- a/src/math/lp/lp_settings.h
+++ b/src/math/lp/lp_settings.h
@@ -339,11 +339,7 @@ struct lp_settings {
m_simplex_strategy = s;
}
-
- bool use_tableau() const {
- return true;
- }
-
+
bool use_tableau_rows() const {
return m_simplex_strategy == simplex_strategy_enum::tableau_rows;
}
diff --git a/src/math/lp/lu.cpp b/src/math/lp/lu.cpp
index 6c9bcc5f65f..313fa990151 100644
--- a/src/math/lp/lu.cpp
+++ b/src/math/lp/lu.cpp
@@ -28,13 +28,10 @@ template double dot_product<double, double>(vector<double> const&, vector<double
template lu<static_matrix<double, double>>::lu(static_matrix<double, double> const&, vector<unsigned int>&, lp_settings&);
template void lu<static_matrix<double, double>>::push_matrix_to_tail(tail_matrix<double, double>*);
template void lu<static_matrix<double, double>>::replace_column(double, indexed_vector<double>&, unsigned);
-template void lu<static_matrix<double, double>>::solve_Bd(unsigned int, indexed_vector<double>&, indexed_vector<double>&);
template lu<static_matrix<double, double>>::~lu();
template void lu<static_matrix<mpq, mpq>>::push_matrix_to_tail(tail_matrix<mpq, mpq>*);
-template void lu<static_matrix<mpq, mpq>>::solve_Bd(unsigned int, indexed_vector<mpq>&, indexed_vector<mpq>&);
template lu<static_matrix<mpq, mpq>>::~lu();
template void lu<static_matrix<mpq, impq>>::push_matrix_to_tail(tail_matrix<mpq, impq >*);
-template void lu<static_matrix<mpq, impq>>::solve_Bd(unsigned int, indexed_vector<mpq>&, indexed_vector<mpq>&);
template lu<static_matrix<mpq, impq>>::~lu();
template mpq dot_product<mpq, mpq>(vector<mpq > const&, vector<mpq > const&);
template void init_factorization<static_matrix<double, double>>
diff --git a/src/math/lp/lu.h b/src/math/lp/lu.h
index aca59065d4b..191018100f8 100644
--- a/src/math/lp/lu.h
+++ b/src/math/lp/lu.h
@@ -171,8 +171,6 @@ class lu {
void print_matrix_compact(std::ostream & f);
void print(indexed_vector<T> & w, const vector<unsigned>& basis);
- void solve_Bd(unsigned a_column, vector<T> & d, indexed_vector<T> & w);
- void solve_Bd(unsigned a_column, indexed_vector<T> & d, indexed_vector<T> & w);
void solve_Bd_faster(unsigned a_column, indexed_vector<T> & d); // d is the right side on the input and the solution at the exit
void solve_yB(vector<T>& y);
diff --git a/src/math/lp/lu_def.h b/src/math/lp/lu_def.h
index 80c9cdf0ed4..059a430124c 100644
--- a/src/math/lp/lu_def.h
+++ b/src/math/lp/lu_def.h
@@ -117,7 +117,7 @@ lu<M>::lu(const M& A,
m_failure(false),
m_row_eta_work_vector(A.row_count()),
m_refactor_counter(0) {
- lp_assert(!(numeric_traits<T>::precise() && settings.use_tableau()));
+ lp_assert(!(numeric_traits<T>::precise() ));
#ifdef Z3DEBUG
debug_test_of_basis(A, basis);
#endif
@@ -256,19 +256,6 @@ void lu< M>::print(indexed_vector<T> & w, const vector<unsigned>& basis) {
print_indexed_vector(w, f);
f.close();
}
-template <typename M>
-void lu< M>::solve_Bd(unsigned a_column, indexed_vector<T> & d, indexed_vector<T> & w) {
- init_vector_w(a_column, w);
-
- if (w.m_index.size() * ratio_of_index_size_to_all_size<T>() < d.m_data.size()) { // this const might need some tuning
- d = w;
- solve_By_for_T_indexed_only(d, m_settings);
- } else {
- d.m_data = w.m_data;
- d.m_index.clear();
- solve_By_when_y_is_ready_for_T(d.m_data, d.m_index);
- }
-}
template <typename M>
void lu< M>::solve_Bd_faster(unsigned a_column, indexed_vector<T> & d) { // puts the a_column into d
diff --git a/src/test/lp/lp.cpp b/src/test/lp/lp.cpp
index fe589835d48..189bd604b92 100644
--- a/src/test/lp/lp.cpp
+++ b/src/test/lp/lp.cpp
@@ -2095,35 +2095,7 @@ void read_indexed_vector(indexed_vector<double> & v, std::ifstream & f) {
}
void check_lu_from_file(std::string lufile_name) {
- std::ifstream f(lufile_name);
- if (!f.is_open()) {
- std::cout << "cannot open file " << lufile_name << std::endl;
- }
- unsigned m, n;
- get_matrix_dimensions(f, m, n);
- std::cout << "init matrix " << m << " by " << n << std::endl;
- static_matrix<double, double> A(m, n);
- read_rows(A, f);
- vector<unsigned> basis;
- read_basis(basis, f);
- indexed_vector<double> v(m);
- // read_indexed_vector(v, f);
- f.close();
- vector<int> basis_heading;
- lp_settings settings;
- vector<unsigned> non_basic_columns;
- lu<static_matrix<double, double>> lsuhl(A, basis, settings);
- indexed_vector<double> d(A.row_count());
- unsigned entering = 26;
- lsuhl.solve_Bd(entering, d, v);
-#ifdef Z3DEBUG
- auto B = get_B(lsuhl, basis);
- vector<double> a(m);
- A.copy_column_to_vector(entering, a);
- indexed_vector<double> cd(d);
- B.apply_from_left(cd.m_data, settings);
- lp_assert(vectors_are_equal(cd.m_data , a));
-#endif
+ lp_assert(false);
}
void test_square_dense_submatrix() {