mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-04-13 12:28:44 +00:00
Code Activity

rm lu

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2023-03-06 07:35:10 -08:00
parent 5f03c93270
commit 62bd3bd1e6
11 changed files with 22 additions and 461 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

241
src/math/lp/lar_core_solver.h
View file

@ -190,16 +190,12 @@ public:
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 @@ public:
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 <typename L, typename K>
void prepare_solver_x_with_signature(const lar_solution_signature & signature, lp_primal_core_solver<L,K> & 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<K>();
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 <typename L, typename K>
void catch_up_in_lu_in_reverse(const vector<unsigned> & trace_of_basis_change, lp_primal_core_solver<L,K> & 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 <typename L, typename K>
void catch_up_in_lu(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading, lp_primal_core_solver<L,K> & 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<L> w(cs.m_A.row_count());
// the queues of delayed indices
std::queue<unsigned> 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<L>(), 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<L>::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<unsigned> & 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 @@ public:
return true;
}
bool catch_up_in_lu_tableau(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading) {
lp_assert(r_basis_is_OK());
// the queues of delayed indices
std::queue<unsigned> 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 @@ public:
}
// returns the trace of basis changes
vector<unsigned> 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<unsigned> 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<unsigned>();
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:

28
src/math/lp/lar_core_solver_def.h
View file

@ -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<unsigned> 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:

11
src/math/lp/lp_core_solver_base.cpp
View file

@ -23,8 +23,6 @@ Revision History:
#include "util/vector.h"
#include <functional>
#include "math/lp/lp_core_solver_base_def.h"
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_when_pivot_row_of_B1_is_ready(unsigned);
template bool lp::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int) const;
@ -33,7 +31,6 @@ template bool lp::lp_core_solver_base<double, double>::find_x_by_solving();
template lp::non_basic_column_value_position lp::lp_core_solver_base<double, double>::get_non_basic_column_value_position(unsigned int) const;
template lp::non_basic_column_value_position lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::get_non_basic_column_value_position(unsigned int) const;
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>&,
vector<unsigned>&,
@ -51,26 +48,20 @@ 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_yB(vector<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_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();
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::find_x_by_solving();
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::init_reduced_costs_for_one_iteration();
template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::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<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_yB(vector<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_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();
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 >&,
@ -106,7 +97,6 @@ template std::string lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq>
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::pretty_print(std::ostream & out);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::restore_state(lp::mpq*, lp::mpq*);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::save_state(lp::mpq*, lp::mpq*);
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::solve_yB(vector<lp::mpq>&) const;
template void lp::lp_core_solver_base<double, double>::init_lu();
template void lp::lp_core_solver_base<lp::mpq, lp::mpq>::init_lu();
template int lp::lp_core_solver_base<double, double>::pivots_in_column_and_row_are_different(int, int) const;
@ -122,7 +112,6 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::column_is_feasible(unsi
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::column_is_feasible(unsigned int) const;
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::snap_non_basic_x_to_bound();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::init_lu();
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::A_mult_x_is_off_on_index(vector<unsigned int> const&) const;
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::find_x_by_solving();
template void lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq> >::restore_x(unsigned int, lp::numeric_pair<lp::mpq> const&);
template bool lp::lp_core_solver_base<lp::mpq, lp::numeric_pair<lp::mpq>>::pivot_column_tableau(unsigned int, unsigned int);

13
src/math/lp/lp_core_solver_base.h
View file

@ -154,9 +154,6 @@ public:
void fill_cb(vector<T> & y) const;
void solve_yB(vector<T> & y) const;
void pretty_print(std::ostream & out);
void save_state(T * w_buffer, T * d_buffer);
@ -175,10 +172,6 @@ public:
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<unsigned> & index) const;
void calculate_pivot_row_when_pivot_row_of_B1_is_ready(unsigned pivot_row);
@ -317,8 +310,6 @@ public:
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 @@ public:
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();

97
src/math/lp/lp_core_solver_base_def.h
View file

@ -116,11 +116,6 @@ fill_cb(vector<T> & y) const {
}
}
template <typename T, typename X> void lp_core_solver_base<T, X>::
solve_yB(vector<T> & 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 <typename T, typename X> void lp_core_solver_base<T, X>::
// update_index_of_ed() {
@ -188,71 +183,6 @@ restore_m_ed(T * buffer) {
}
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
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);
if (delta != numeric_traits<X>::zero()) {
return true;
}
}
return false;
}
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));
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 <typename T, typename X> bool lp_core_solver_base<T, X>::
A_mult_x_is_off_on_index(const vector<unsigned> & index) const {
lp_assert(m_x.size() == m_A.column_count());
if (numeric_traits<T>::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<X>::zero()) {
return true;
}
}
return false;
#endif
// todo(levnach) run on m_ed.m_index only !!!!!
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));
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 <typename T, typename X> void lp_core_solver_base<T, X>::
@ -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 <typename T, typename X> bool lp_core_solver_base<T, X>::
@ -411,7 +341,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();
return !A_mult_x_is_off();
return true;
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::column_is_feasible(unsigned j) const {
@ -591,24 +521,6 @@ template <typename T, typename X> bool lp_core_solver_base<T, X>::
return true;
}
template <typename T, typename X> void lp_core_solver_base<T, X>::
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 <typename T, typename X> void lp_core_solver_base<T, X>::
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 <typename T, typename X> void lp_core_solver_base<T, X>::
init_reduced_costs_for_one_iteration() {
solve_yB(m_y);
fill_reduced_costs_from_m_y_by_rows();
}
template <typename T, typename X> non_basic_column_value_position lp_core_solver_base<T, X>::
get_non_basic_column_value_position(unsigned j) const {

2
src/math/lp/lp_dual_core_solver.h
View file

@ -76,7 +76,7 @@ public:
m_a_wave(this->m_m()),
m_betas(this->m_m()) {
m_harris_tolerance = numeric_traits<T>::precise()? numeric_traits<T>::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();
}

14
src/math/lp/lp_dual_core_solver_def.h
View file

@ -74,8 +74,7 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::recalculate_xB
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::recalculate_d() {
this->solve_yB(this->m_y);
this->fill_reduced_costs_from_m_y_by_rows();
lp_assert(false)
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::init_betas() {
@ -316,15 +315,8 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_d() {
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::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<T>::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;
}

5
src/math/lp/lp_primal_core_solver_def.h
View file

@ -705,10 +705,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
return 0;
}
if ((!numeric_traits<T>::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();

9
src/math/lp/lp_primal_core_solver_tableau_def.h
View file

@ -107,10 +107,6 @@ unsigned lp_primal_core_solver<T, X>::solve_with_tableau() {
return 0;
}
if ((!numeric_traits<T>::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<T, X>::solve_with_tableau() {
}
template <typename T, typename X>void lp_primal_core_solver<T, X>::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 <typename T, typename X>void lp_primal_core_solver<T, X>::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 <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
template <typename T, typename X>
void lp_primal_core_solver<T, X>::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 <typename T, typename X> int lp_primal_core_solver<T, X>::find_leaving_
return m_leaving_candidates[k];
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::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 <typename T, typename X> void lp_primal_core_solver<T, X>::

3
src/math/lp/lu.h
View file

@ -173,8 +173,7 @@ public:
void print(indexed_vector<T> & w, const vector<unsigned>& basis);
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);
void solve_yB_indexed(indexed_vector<T>& y);
void add_delta_to_solution_indexed(indexed_vector<T>& y);

60
src/math/lp/lu_def.h
View file

@ -263,20 +263,6 @@ void lu< M>::solve_Bd_faster(unsigned a_column, indexed_vector<T> & d) { // puts
solve_By_for_T_indexed_only(d, m_settings);
}
template <typename M>
void lu< M>::solve_yB(vector<T>& 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 <typename M>
void lu< M>::solve_yB_indexed(indexed_vector<T>& y) {
lp_assert(y.is_OK());
@ -412,55 +398,15 @@ void lu< M>::find_error_of_yB_indexed(const indexed_vector<T>& y, const vector<i
// y is the input
template <typename M>
void lu< M>::solve_yB_with_error_check_indexed(indexed_vector<T> & y, const vector<int>& heading, const vector<unsigned> & basis, const lp_settings & settings) {
if (numeric_traits<T>::precise()) {
if (y.m_index.size() * ratio_of_index_size_to_all_size<T>() * 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<T>() < 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<T>() >= 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 <typename M>
void lu< M>::solve_yB_with_error_check(vector<T> & y, const vector<unsigned>& basis) {
if (numeric_traits<T>::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 <typename M>
void lu< M>::apply_Q_R_to_U(permutation_matrix<T, X> & r_wave) {