diff --git a/src/math/lp/lar_core_solver.h b/src/math/lp/lar_core_solver.h index 9678edd6bed..999eef13bc8 100644 --- a/src/math/lp/lar_core_solver.h +++ b/src/math/lp/lar_core_solver.h @@ -190,16 +190,12 @@ class lar_core_solver { m_r_upper_bounds.pop(k); m_column_types.pop(k); - delete m_r_solver.m_factorization; - m_r_solver.m_factorization = nullptr; 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()); m_d_A.pop(k); // doubles - delete m_d_solver.m_factorization; - m_d_solver.m_factorization = nullptr; m_d_x.resize(m_d_A.column_count()); pop_basis(k); @@ -294,172 +290,13 @@ class lar_core_solver { unsigned jb = m_r_solver.m_basis[i]; m_r_solver.add_delta_to_x_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc)); } - CASSERT("A_off", m_r_solver.A_mult_x_is_off() == false); + } lp_assert(m_r_solver.inf_set_is_correct()); } - template - void prepare_solver_x_with_signature(const lar_solution_signature & signature, lp_primal_core_solver & s) { - for (auto &t : signature) { - unsigned j = t.first; - lp_assert(m_r_heading[j] < 0); - auto pos_type = t.second; - switch (pos_type) { - case at_lower_bound: - s.m_x[j] = s.m_lower_bounds[j]; - break; - case at_fixed: - case at_upper_bound: - s.m_x[j] = s.m_upper_bounds[j]; - break; - case free_of_bounds: { - s.m_x[j] = zero_of_type(); - continue; - } - case not_at_bound: - 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; - case column_type::lower_bound: - s.m_x[j] = s.m_lower_bounds[j]; - break; - case column_type::boxed: - if (settings().random_next() % 2) { - s.m_x[j] = s.m_lower_bounds[j]; - } else { - s.m_x[j] = s.m_upper_bounds[j]; - } - break; - case column_type::fixed: - s.m_x[j] = s.m_lower_bounds[j]; - break; - default: - lp_assert(false); - } - break; - default: - lp_unreachable(); - } - } - - // lp_assert(is_zero_vector(s.m_b)); - s.solve_Ax_eq_b(); - } - - template - void catch_up_in_lu_in_reverse(const vector & trace_of_basis_change, lp_primal_core_solver & cs) { - // recover the previous working basis - for (unsigned i = trace_of_basis_change.size(); i > 0; i-= 2) { - unsigned entering = trace_of_basis_change[i-1]; - unsigned leaving = trace_of_basis_change[i-2]; - cs.change_basis_unconditionally(entering, leaving); - } - cs.init_lu(); - } - - //basis_heading is the basis heading of the solver owning trace_of_basis_change - // here we compact the trace as we go to avoid unnecessary column changes - template - void catch_up_in_lu(const vector & trace_of_basis_change, const vector & basis_heading, lp_primal_core_solver & cs) { - if (cs.m_factorization == nullptr || cs.m_factorization->m_refactor_counter + trace_of_basis_change.size()/2 >= 200) { - for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { - unsigned entering = trace_of_basis_change[i]; - unsigned leaving = trace_of_basis_change[i+1]; - cs.change_basis_unconditionally(entering, leaving); - } - if (cs.m_factorization != nullptr) { - delete cs.m_factorization; - cs.m_factorization = nullptr; - } - } else { - indexed_vector w(cs.m_A.row_count()); - // the queues of delayed indices - std::queue entr_q, leav_q; - auto * l = cs.m_factorization; - lp_assert(l->get_status() == LU_status::OK); - for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { - unsigned entering = trace_of_basis_change[i]; - unsigned leaving = trace_of_basis_change[i+1]; - bool good_e = basis_heading[entering] >= 0 && cs.m_basis_heading[entering] < 0; - bool good_l = basis_heading[leaving] < 0 && cs.m_basis_heading[leaving] >= 0; - if (!good_e && !good_l) continue; - if (good_e && !good_l) { - while (!leav_q.empty() && cs.m_basis_heading[leav_q.front()] < 0) - leav_q.pop(); - if (!leav_q.empty()) { - leaving = leav_q.front(); - leav_q.pop(); - } else { - entr_q.push(entering); - continue; - } - } else if (!good_e && good_l) { - while (!entr_q.empty() && cs.m_basis_heading[entr_q.front()] >= 0) - entr_q.pop(); - if (!entr_q.empty()) { - entering = entr_q.front(); - entr_q.pop(); - } else { - leav_q.push(leaving); - continue; - } - } - lp_assert(cs.m_basis_heading[entering] < 0); - lp_assert(cs.m_basis_heading[leaving] >= 0); - if (l->get_status() == LU_status::OK) { - l->prepare_entering(entering, w); // to init vector w - l->replace_column(zero_of_type(), w, cs.m_basis_heading[leaving]); - } - cs.change_basis_unconditionally(entering, leaving); - } - if (l->get_status() != LU_status::OK) { - delete l; - cs.m_factorization = nullptr; - } - } - if (cs.m_factorization == nullptr) { - if (numeric_traits::precise()) - init_factorization(cs.m_factorization, cs.m_A, cs.m_basis, settings()); - } - } - - bool no_r_lu() const { - return m_r_solver.m_factorization == nullptr || m_r_solver.m_factorization->get_status() == LU_status::Degenerated; - } - - void solve_on_signature_tableau(const lar_solution_signature & signature, const vector & changes_of_basis) { - r_basis_is_OK(); - 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 - prepare_solver_x_with_signature_tableau(signature); // still are going to use the signature partially - 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_tableau(signature); - 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); - } - lp_assert(r_basis_is_OK()); - } - + bool adjust_x_of_column(unsigned j) { /* if (m_r_solver.m_basis_heading[j] >= 0) { @@ -478,58 +315,6 @@ class lar_core_solver { return true; } - - bool catch_up_in_lu_tableau(const vector & trace_of_basis_change, const vector & basis_heading) { - lp_assert(r_basis_is_OK()); - // the queues of delayed indices - std::queue entr_q, leav_q; - for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { - unsigned entering = trace_of_basis_change[i]; - unsigned leaving = trace_of_basis_change[i+1]; - bool good_e = basis_heading[entering] >= 0 && m_r_solver.m_basis_heading[entering] < 0; - bool good_l = basis_heading[leaving] < 0 && m_r_solver.m_basis_heading[leaving] >= 0; - if (!good_e && !good_l) continue; - if (good_e && !good_l) { - while (!leav_q.empty() && m_r_solver.m_basis_heading[leav_q.front()] < 0) - leav_q.pop(); - if (!leav_q.empty()) { - leaving = leav_q.front(); - leav_q.pop(); - } else { - entr_q.push(entering); - continue; - } - } else if (!good_e && good_l) { - while (!entr_q.empty() && m_r_solver.m_basis_heading[entr_q.front()] >= 0) - entr_q.pop(); - if (!entr_q.empty()) { - entering = entr_q.front(); - entr_q.pop(); - } else { - leav_q.push(leaving); - continue; - } - } - lp_assert(m_r_solver.m_basis_heading[entering] < 0); - lp_assert(m_r_solver.m_basis_heading[leaving] >= 0); - m_r_solver.change_basis_unconditionally(entering, leaving); - if(!m_r_solver.pivot_column_tableau(entering, m_r_solver.m_basis_heading[entering])) { - // unroll the last step - m_r_solver.change_basis_unconditionally(leaving, entering); -#ifdef Z3DEBUG - bool t = -#endif - m_r_solver.pivot_column_tableau(leaving, m_r_solver.m_basis_heading[leaving]); -#ifdef Z3DEBUG - lp_assert(t); -#endif - return false; - } - } - lp_assert(r_basis_is_OK()); - return true; - } - bool r_basis_is_OK() const { #ifdef Z3DEBUG @@ -598,27 +383,7 @@ class lar_core_solver { } - // returns the trace of basis changes - vector find_solution_signature_with_doubles(lar_solution_signature & signature) { - if (m_d_solver.m_factorization == nullptr || m_d_solver.m_factorization->get_status() != LU_status::OK) { - vector ret; - return ret; - } - get_bounds_for_double_solver(); - - extract_signature_from_lp_core_solver(m_r_solver, signature); - prepare_solver_x_with_signature(signature, m_d_solver); - m_d_solver.start_tracing_basis_changes(); - m_d_solver.find_feasible_solution(); - if (settings().get_cancel_flag()) - return vector(); - - m_d_solver.stop_tracing_basis_changes(); - extract_signature_from_lp_core_solver(m_d_solver, signature); - return m_d_solver.m_trace_of_basis_change_vector; - } - - + bool lower_bound_is_set(unsigned j) const { switch (m_column_types[j]) { case column_type::free_column: diff --git a/src/math/lp/lar_core_solver_def.h b/src/math/lp/lar_core_solver_def.h index 182029e3e98..22dc23bd512 100644 --- a/src/math/lp/lar_core_solver_def.h +++ b/src/math/lp/lar_core_solver_def.h @@ -78,7 +78,6 @@ void lar_core_solver::prefix_d() { } void lar_core_solver::fill_not_improvable_zero_sum_from_inf_row() { - CASSERT("A_off", m_r_solver.A_mult_x_is_off() == false); 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(); @@ -127,29 +126,16 @@ void lar_core_solver::solve() { return; } ++settings().stats().m_need_to_solve_inf; - CASSERT("A_off", !m_r_solver.A_mult_x_is_off()); lp_assert( r_basis_is_OK()); - if (need_to_presolve_with_double_solver()) { - TRACE("lar_solver", tout << "presolving\n";); - prefix_d(); - lar_solution_signature solution_signature; - vector changes_of_basis = find_solution_signature_with_doubles(solution_signature); - if (m_d_solver.get_status() == lp_status::TIME_EXHAUSTED) { - m_r_solver.set_status(lp_status::TIME_EXHAUSTED); - return; - } - solve_on_signature_tableau(solution_signature, changes_of_basis); - - lp_assert( r_basis_is_OK()); - } else { + - 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(r_basis_is_OK()); + 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(r_basis_is_OK()); + switch (m_r_solver.get_status()) { case lp_status::INFEASIBLE: diff --git a/src/math/lp/lp_core_solver_base.cpp b/src/math/lp/lp_core_solver_base.cpp index 9f6f2534f6c..508a06d730e 100644 --- a/src/math/lp/lp_core_solver_base.cpp +++ b/src/math/lp/lp_core_solver_base.cpp @@ -23,8 +23,6 @@ Revision History: #include "util/vector.h" #include #include "math/lp/lp_core_solver_base_def.h" -template bool lp::lp_core_solver_base::A_mult_x_is_off() const; -template bool lp::lp_core_solver_base::A_mult_x_is_off_on_index(const vector &) const; template bool lp::lp_core_solver_base::basis_heading_is_correct() const; template void lp::lp_core_solver_base::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned); template bool lp::lp_core_solver_base::column_is_dual_feasible(unsigned int) const; @@ -33,7 +31,6 @@ template bool lp::lp_core_solver_base::find_x_by_solving(); template lp::non_basic_column_value_position lp::lp_core_solver_base::get_non_basic_column_value_position(unsigned int) const; template lp::non_basic_column_value_position lp::lp_core_solver_base >::get_non_basic_column_value_position(unsigned int) const; template lp::non_basic_column_value_position lp::lp_core_solver_base::get_non_basic_column_value_position(unsigned int) const; -template void lp::lp_core_solver_base::init_reduced_costs_for_one_iteration(); template lp::lp_core_solver_base::lp_core_solver_base( lp::static_matrix&, // vector&, vector&, @@ -51,26 +48,20 @@ template void lp::lp_core_solver_base::set_non_basic_x_to_correc template void lp::lp_core_solver_base::snap_xN_to_bounds_and_free_columns_to_zeroes(); template void lp::lp_core_solver_base >::snap_xN_to_bounds_and_free_columns_to_zeroes(); template void lp::lp_core_solver_base::solve_Ax_eq_b(); -template void lp::lp_core_solver_base::solve_yB(vector&) const; template void lp::lp_core_solver_base::add_delta_to_entering(unsigned int, const double&); -template bool lp::lp_core_solver_base::A_mult_x_is_off() const; -template bool lp::lp_core_solver_base::A_mult_x_is_off_on_index(const vector &) const; template bool lp::lp_core_solver_base::basis_heading_is_correct() const ; template void lp::lp_core_solver_base::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned); template bool lp::lp_core_solver_base::column_is_dual_feasible(unsigned int) const; template void lp::lp_core_solver_base::fill_reduced_costs_from_m_y_by_rows(); template bool lp::lp_core_solver_base::find_x_by_solving(); -template void lp::lp_core_solver_base::init_reduced_costs_for_one_iteration(); template bool lp::lp_core_solver_base::print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over(char const*, std::ostream &); template void lp::lp_core_solver_base::restore_x(unsigned int, lp::mpq const&); template void lp::lp_core_solver_base::set_non_basic_x_to_correct_bounds(); template void lp::lp_core_solver_base::solve_Ax_eq_b(); -template void lp::lp_core_solver_base::solve_yB(vector&) const; template void lp::lp_core_solver_base::add_delta_to_entering(unsigned int, const lp::mpq&); template void lp::lp_core_solver_base >::calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned); template void lp::lp_core_solver_base >::init(); template void lp::lp_core_solver_base >::init_basis_heading_and_non_basic_columns_vector(); -template void lp::lp_core_solver_base >::init_reduced_costs_for_one_iteration(); template lp::lp_core_solver_base >::lp_core_solver_base(lp::static_matrix >&, // vector >&, vector&, vector &, vector &, vector >&, vector&, lp::lp_settings&, const column_namer&, const vector&, @@ -106,7 +97,6 @@ template std::string lp::lp_core_solver_base template void lp::lp_core_solver_base >::pretty_print(std::ostream & out); template void lp::lp_core_solver_base >::restore_state(lp::mpq*, lp::mpq*); template void lp::lp_core_solver_base >::save_state(lp::mpq*, lp::mpq*); -template void lp::lp_core_solver_base >::solve_yB(vector&) const; template void lp::lp_core_solver_base::init_lu(); template void lp::lp_core_solver_base::init_lu(); template int lp::lp_core_solver_base::pivots_in_column_and_row_are_different(int, int) const; @@ -122,7 +112,6 @@ template bool lp::lp_core_solver_base::column_is_feasible(unsi template bool lp::lp_core_solver_base >::column_is_feasible(unsigned int) const; template bool lp::lp_core_solver_base >::snap_non_basic_x_to_bound(); template void lp::lp_core_solver_base >::init_lu(); -template bool lp::lp_core_solver_base >::A_mult_x_is_off_on_index(vector const&) const; template bool lp::lp_core_solver_base >::find_x_by_solving(); template void lp::lp_core_solver_base >::restore_x(unsigned int, lp::numeric_pair const&); template bool lp::lp_core_solver_base>::pivot_column_tableau(unsigned int, unsigned int); diff --git a/src/math/lp/lp_core_solver_base.h b/src/math/lp/lp_core_solver_base.h index dce805b62ee..edff349d23d 100644 --- a/src/math/lp/lp_core_solver_base.h +++ b/src/math/lp/lp_core_solver_base.h @@ -154,9 +154,6 @@ class lp_core_solver_base { void fill_cb(vector & y) const; - void solve_yB(vector & y) const; - - void pretty_print(std::ostream & out); void save_state(T * w_buffer, T * d_buffer); @@ -175,10 +172,6 @@ class lp_core_solver_base { void copy_m_ed(T * buffer); void restore_m_ed(T * buffer); - - bool A_mult_x_is_off() const; - - bool A_mult_x_is_off_on_index(const vector & index) const; void calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned pivot_row); @@ -317,8 +310,6 @@ class lp_core_solver_base { bool basis_heading_is_correct() const; - void restore_x_and_refactor(int entering, int leaving, X const & t); - void restore_x(unsigned entering, X const & t); void fill_reduced_costs_from_m_y_by_rows(); @@ -435,9 +426,7 @@ class lp_core_solver_base { void snap_xN_to_bounds_and_fill_xB(); void snap_xN_to_bounds_and_free_columns_to_zeroes(); - - void init_reduced_costs_for_one_iteration(); - + non_basic_column_value_position get_non_basic_column_value_position(unsigned j) const; void init_lu(); diff --git a/src/math/lp/lp_core_solver_base_def.h b/src/math/lp/lp_core_solver_base_def.h index 51b24128fef..ae04bd5eec6 100644 --- a/src/math/lp/lp_core_solver_base_def.h +++ b/src/math/lp/lp_core_solver_base_def.h @@ -116,11 +116,6 @@ fill_cb(vector & y) const { } } -template void lp_core_solver_base:: -solve_yB(vector & y) const { - fill_cb(y); // now y = cB, that is the projection of costs to basis - m_factorization->solve_yB_with_error_check(y, m_basis); -} // template void lp_core_solver_base:: // update_index_of_ed() { @@ -188,71 +183,6 @@ restore_m_ed(T * buffer) { } } -template bool lp_core_solver_base:: -A_mult_x_is_off() const { - lp_assert(m_x.size() == m_A.column_count()); - if (numeric_traits::precise()) { - for (unsigned i = 0; i < m_m(); i++) { - X delta = /*m_b[i] */- m_A.dot_product_with_row(i, m_x); - if (delta != numeric_traits::zero()) { - return true; - } - } - return false; - } - T feps = convert_struct::convert(m_settings.refactor_tolerance); - X one = convert_struct::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)); - auto eps = feps /* * (one + T(0.1) * abs(m_b[i])) */; - - if (delta > eps) { -#if 0 - LP_OUT(m_settings, "x is off (" - << "m_b[" << i << "] = " << m_b[i] << " " - << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ' - << "delta = " << delta << ' ' - << "iters = " << total_iterations() << ")" << std::endl); -#endif - return true; - } - } - return false; -} -template bool lp_core_solver_base:: -A_mult_x_is_off_on_index(const vector & index) const { - lp_assert(m_x.size() == m_A.column_count()); - if (numeric_traits::precise()) return false; -#if RUN_A_MULT_X_IS_OFF_FOR_PRECESE - for (unsigned i : index) { - X delta = m_b[i] - m_A.dot_product_with_row(i, m_x); - if (delta != numeric_traits::zero()) { - return true; - } - } - return false; -#endif - // todo(levnach) run on m_ed.m_index only !!!!! - T feps = convert_struct::convert(m_settings.refactor_tolerance); - X one = convert_struct::convert(1.0); - for (unsigned i : index) { - X delta = abs(/*m_b[i] -*/ m_A.dot_product_with_row(i, m_x)); - auto eps = feps /* *(one + T(0.1) * abs(m_b[i])) */; - - if (delta > eps) { -#if 0 - LP_OUT(m_settings, "x is off (" - << "m_b[" << i << "] = " << m_b[i] << " " - << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ' - << "delta = " << delta << ' ' - << "iters = " << total_iterations() << ")" << std::endl); -#endif - return true; - } - } - return false; -} - template void lp_core_solver_base:: @@ -292,7 +222,7 @@ print_statistics(char const* str, X cost, std::ostream & out) { if (str!= nullptr) out << str << " "; out << "iterations = " << (total_iterations() - 1) << ", cost = " << T_to_string(cost) - << ", nonzeros = " << (m_factorization != nullptr? m_factorization->get_number_of_nonzeroes() : m_A.number_of_non_zeroes()) << std::endl; + << ", nonzeros = " << m_A.number_of_non_zeroes() << std::endl; } template bool lp_core_solver_base:: @@ -411,7 +341,7 @@ rs_minus_Anx(vector & rs) { template bool lp_core_solver_base:: find_x_by_solving() { solve_Ax_eq_b(); - return !A_mult_x_is_off(); + return true; } template bool lp_core_solver_base::column_is_feasible(unsigned j) const { @@ -591,24 +521,6 @@ template bool lp_core_solver_base:: return true; } -template void lp_core_solver_base:: -restore_x_and_refactor(int entering, int leaving, X const & t) { - this->restore_basis_change(entering, leaving); - restore_x(entering, t); - init_factorization(m_factorization, m_A, m_basis, m_settings); - if (m_factorization->get_status() == LU_status::Degenerated) { - LP_OUT(m_settings, "cannot refactor" << std::endl); - m_status = lp_status::FLOATING_POINT_ERROR; - return; - } - // solve_Ax_eq_b(); - if (A_mult_x_is_off()) { - LP_OUT(m_settings, "cannot restore solution" << std::endl); - m_status = lp_status::FLOATING_POINT_ERROR; - return; - } -} - template void lp_core_solver_base:: restore_x(unsigned entering, X const & t) { if (is_zero(t)) return; @@ -731,11 +643,6 @@ snap_xN_to_bounds_and_free_columns_to_zeroes() { solve_Ax_eq_b(); } -template void lp_core_solver_base:: -init_reduced_costs_for_one_iteration() { - solve_yB(m_y); - fill_reduced_costs_from_m_y_by_rows(); -} template non_basic_column_value_position lp_core_solver_base:: get_non_basic_column_value_position(unsigned j) const { diff --git a/src/math/lp/lp_dual_core_solver.h b/src/math/lp/lp_dual_core_solver.h index 804879c3da2..b78a2214701 100644 --- a/src/math/lp/lp_dual_core_solver.h +++ b/src/math/lp/lp_dual_core_solver.h @@ -76,7 +76,7 @@ class lp_dual_core_solver:public lp_core_solver_base { m_a_wave(this->m_m()), m_betas(this->m_m()) { m_harris_tolerance = numeric_traits::precise()? numeric_traits::zero() : T(this->m_settings.harris_feasibility_tolerance); - this->solve_yB(this->m_y); + lp_assert(false); this->init_basic_part_of_basis_heading(); fill_non_basis_with_only_able_to_enter_columns(); } diff --git a/src/math/lp/lp_dual_core_solver_def.h b/src/math/lp/lp_dual_core_solver_def.h index fd8fc80718f..9c39fe195bf 100644 --- a/src/math/lp/lp_dual_core_solver_def.h +++ b/src/math/lp/lp_dual_core_solver_def.h @@ -74,8 +74,7 @@ template void lp_dual_core_solver::recalculate_xB } template void lp_dual_core_solver::recalculate_d() { - this->solve_yB(this->m_y); - this->fill_reduced_costs_from_m_y_by_rows(); +lp_assert(false) } template void lp_dual_core_solver::init_betas() { @@ -316,15 +315,8 @@ template void lp_dual_core_solver::restore_d() { } template bool lp_dual_core_solver::d_is_correct() { - this->solve_yB(this->m_y); - for (auto j : this->non_basis()) { - T d = this->m_costs[j] - this->m_A.dot_product_with_column(this->m_y, j); - if (numeric_traits::get_double(abs(d - this->m_d[j])) >= 0.001) { - LP_OUT(this->m_settings, "total_iterations = " << this->total_iterations() << std::endl - << "d[" << j << "] = " << this->m_d[j] << " but should be " << d << std::endl); - return false; - } - } + + lp_assert(false); 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 7d522f9331d..19ffaa1eadc 100644 --- a/src/math/lp/lp_primal_core_solver_def.h +++ b/src/math/lp/lp_primal_core_solver_def.h @@ -705,10 +705,7 @@ template unsigned lp_primal_core_solver::solve() return 0; } - if ((!numeric_traits::precise()) && this->A_mult_x_is_off()) { - this->set_status(lp_status::FLOATING_POINT_ERROR); - return 0; - } + do { if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->using_infeas_costs()? "inf" : "feas"), * this->m_settings.get_message_ostream())) { return this->total_iterations(); 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 fa25694adc7..736b79db5ad 100644 --- a/src/math/lp/lp_primal_core_solver_tableau_def.h +++ b/src/math/lp/lp_primal_core_solver_tableau_def.h @@ -107,10 +107,6 @@ unsigned lp_primal_core_solver::solve_with_tableau() { return 0; } - if ((!numeric_traits::precise()) && this->A_mult_x_is_off()) { - this->set_status(lp_status::FLOATING_POINT_ERROR); - return 0; - } do { if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over((this->using_infeas_costs()? "inf t" : "feas t"), * this->m_settings.get_message_ostream())) { return this->total_iterations(); @@ -183,7 +179,6 @@ unsigned lp_primal_core_solver::solve_with_tableau() { } template void lp_primal_core_solver::advance_on_entering_and_leaving_tableau(int entering, int leaving, X & t) { - CASSERT("A_off", this->A_mult_x_is_off() == false); lp_assert(leaving >= 0 && entering >= 0); lp_assert((this->m_settings.simplex_strategy() == simplex_strategy_enum::tableau_rows) || @@ -200,7 +195,6 @@ template void lp_primal_core_solver::advance_on_en t = -t; } this->update_basis_and_x_tableau(entering, leaving, t); - CASSERT("A_off", this->A_mult_x_is_off() == false); this->iters_with_no_cost_growing() = 0; } else { this->pivot_column_tableau(entering, this->m_basis_heading[leaving]); @@ -224,7 +218,6 @@ template void lp_primal_core_solver::advance_on_en template void lp_primal_core_solver::advance_on_entering_equal_leaving_tableau(int entering, X & t) { - CASSERT("A_off", !this->A_mult_x_is_off() ); this->update_x_tableau(entering, t * m_sign_of_entering_delta); if (this->m_look_for_feasible_solution_only && this->current_x_is_feasible()) return; @@ -296,7 +289,6 @@ template int lp_primal_core_solver::find_leaving_ return m_leaving_candidates[k]; } template void lp_primal_core_solver::init_run_tableau() { - CASSERT("A_off", this->A_mult_x_is_off() == false); lp_assert(basis_columns_are_set_correctly()); this->m_basis_sort_counter = 0; // to initiate the sort of the basis // this->set_total_iterations(0); @@ -350,7 +342,6 @@ update_x_tableau(unsigned entering, const X& delta) { this->insert_column_into_inf_set(j); } } - CASSERT("A_off", this->A_mult_x_is_off() == false); } template void lp_primal_core_solver:: diff --git a/src/math/lp/lu.h b/src/math/lp/lu.h index b3e2f0c2221..18251e3ada5 100644 --- a/src/math/lp/lu.h +++ b/src/math/lp/lu.h @@ -173,8 +173,7 @@ class lu { void print(indexed_vector & w, const vector& basis); void solve_Bd_faster(unsigned a_column, indexed_vector & d); // d is the right side on the input and the solution at the exit - void solve_yB(vector& y); - + void solve_yB_indexed(indexed_vector& y); void add_delta_to_solution_indexed(indexed_vector& y); diff --git a/src/math/lp/lu_def.h b/src/math/lp/lu_def.h index d7c0c0c3a07..ca34481cca8 100644 --- a/src/math/lp/lu_def.h +++ b/src/math/lp/lu_def.h @@ -263,20 +263,6 @@ void lu< M>::solve_Bd_faster(unsigned a_column, indexed_vector & d) { // puts solve_By_for_T_indexed_only(d, m_settings); } -template -void lu< M>::solve_yB(vector& y) { - // first solve yU = cb*R(-1) - m_R.apply_reverse_from_right_to_T(y); // got y = cb*R(-1) - m_U.solve_y_U(y); // got y*U=cb*R(-1) - m_Q.apply_reverse_from_right_to_T(y); // - for (auto e = m_tail.rbegin(); e != m_tail.rend(); ++e) { -#ifdef Z3DEBUG - (*e)->set_number_of_columns(m_dim); -#endif - (*e)->apply_from_right(y); - } -} - template void lu< M>::solve_yB_indexed(indexed_vector& y) { lp_assert(y.is_OK()); @@ -412,55 +398,15 @@ void lu< M>::find_error_of_yB_indexed(const indexed_vector& y, const vector void lu< M>::solve_yB_with_error_check_indexed(indexed_vector & y, const vector& heading, const vector & basis, const lp_settings & settings) { - if (numeric_traits::precise()) { - if (y.m_index.size() * ratio_of_index_size_to_all_size() * 3 < m_A.column_count()) { - solve_yB_indexed(y); - } else { - solve_yB(y.m_data); - y.restore_index_and_clean_from_data(); - } - return; - } - lp_assert(m_y_copy.is_OK()); - lp_assert(y.is_OK()); - if (y.m_index.size() * ratio_of_index_size_to_all_size() < m_A.column_count()) { - m_y_copy = y; - solve_yB_indexed(y); - lp_assert(y.is_OK()); - if (y.m_index.size() * ratio_of_index_size_to_all_size() >= m_A.column_count()) { - find_error_of_yB(m_y_copy.m_data, y.m_data, basis); - solve_yB(m_y_copy.m_data); - add_delta_to_solution(m_y_copy.m_data, y.m_data); - y.restore_index_and_clean_from_data(); - m_y_copy.clear_all(); - } else { - find_error_of_yB_indexed(y, heading, settings); // this works with m_y_copy - solve_yB_indexed(m_y_copy); - add_delta_to_solution_indexed(y); - } - lp_assert(m_y_copy.is_OK()); - } else { - solve_yB_with_error_check(y.m_data, basis); - y.restore_index_and_clean_from_data(); - } -} + lp_assert(false); + } // solves y*B = y // y is the input template void lu< M>::solve_yB_with_error_check(vector & y, const vector& basis) { - if (numeric_traits::precise()) { - solve_yB(y); - return; - } - auto & yc = m_y_copy.m_data; - yc =y; // copy y aside - solve_yB(y); - find_error_of_yB(yc, y, basis); - solve_yB(yc); - add_delta_to_solution(yc, y); - m_y_copy.clear_all(); + lp_assert(false); } template void lu< M>::apply_Q_R_to_U(permutation_matrix & r_wave) {