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_lp

Browse source 
Signed-off-by: Lev Nachmanson <levnach@hotmail.com>
This commit is contained in:
Lev Nachmanson 2023-03-04 14:58:49 -08:00
parent 527f0d1242
commit 25f103db1a
19 changed files with 65 additions and 641 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/math/lp/core_solver_pretty_printer.h
View file

@ -110,7 +110,6 @@ public:
return T_to_string(m_exact_column_norms[col]);
}
void print_exact_norms();
void print_approx_norms();

67
src/math/lp/core_solver_pretty_printer_def.h
View file

@ -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: ";

1
src/math/lp/int_solver.cpp
View file

@ -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";);

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

@ -80,8 +80,7 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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];

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

@ -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())
{

57
src/math/lp/lar_solver.cpp
View file

@ -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();

22
src/math/lp/lar_solver.h
View file

@ -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 @@ public:
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 @@ public:
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 @@ public:
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,

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

@ -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;

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

@ -156,11 +156,6 @@ public:
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 @@ public:
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 @@ public:
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 @@ public:
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];
}

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

@ -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);
}
}

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

@ -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() {

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

@ -468,7 +468,6 @@ public:
}
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 @@ public:
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; }

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

@ -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) {

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

@ -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]);

6
src/math/lp/lp_settings.h
View file

@ -339,11 +339,7 @@ public:
m_simplex_strategy = s;
}
bool use_tableau() const {
return true;
}
bool use_tableau_rows() const {
return m_simplex_strategy == simplex_strategy_enum::tableau_rows;
}

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

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

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

@ -171,8 +171,6 @@ public:
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);

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

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

30
src/test/lp/lp.cpp
View file

@ -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() {