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
This commit is contained in:
Lev Nachmanson 2023-03-06 07:44:22 -08:00
parent 62bd3bd1e6
commit 1da4c018e4
8 changed files with 0 additions and 969 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/math/lp/CMakeLists.txt
View file

@ -19,7 +19,6 @@ z3_add_component(lp
lar_solver.cpp
lar_core_solver.cpp
lp_core_solver_base.cpp
lp_dual_core_solver.cpp
lp_primal_core_solver.cpp
lp_settings.cpp
lu.cpp

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

@ -27,7 +27,6 @@ template bool lp::lp_core_solver_base<double, double>::basis_heading_is_correct(
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();
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;
@ -53,7 +52,6 @@ template bool lp::lp_core_solver_base<lp::mpq, lp::mpq>::basis_heading_is_correc
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 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();
@ -112,7 +110,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> >::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);
template bool lp::lp_core_solver_base<double, double>::pivot_column_tableau(unsigned int, unsigned int);

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

@ -298,9 +298,6 @@ public:
void rs_minus_Anx(vector<X> & rs);
bool find_x_by_solving();
bool basis_has_no_doubles() const;
bool non_basis_has_no_doubles() const;

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

@ -338,12 +338,6 @@ 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 true;
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::column_is_feasible(unsigned j) const {
const X& x = this->m_x[j];
switch (this->m_column_types[j]) {

44
src/math/lp/lp_dual_core_solver.cpp
View file

@ -1,44 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#include <utility>
#include <memory>
#include <string>
#include "util/vector.h"
#include <functional>
#include "math/lp/lp_dual_core_solver_def.h"
template void lp::lp_dual_core_solver<lp::mpq, lp::mpq>::start_with_initial_basis_and_make_it_dual_feasible();
template void lp::lp_dual_core_solver<lp::mpq, lp::mpq>::solve();
template lp::lp_dual_core_solver<double, double>::lp_dual_core_solver(lp::static_matrix<double, double>&, vector<bool>&,
vector<double>&,
vector<double>&,
vector<unsigned int>&,
vector<unsigned> &,
vector<int> &,
vector<double>&,
vector<lp::column_type>&,
vector<double>&,
vector<double>&,
lp::lp_settings&, const lp::column_namer&);
template void lp::lp_dual_core_solver<double, double>::start_with_initial_basis_and_make_it_dual_feasible();
template void lp::lp_dual_core_solver<double, double>::solve();
template void lp::lp_dual_core_solver<lp::mpq, lp::mpq>::restore_non_basis();
template void lp::lp_dual_core_solver<double, double>::restore_non_basis();
template void lp::lp_dual_core_solver<double, double>::revert_to_previous_basis();
template void lp::lp_dual_core_solver<lp::mpq, lp::mpq>::revert_to_previous_basis();

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

@ -1,212 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include "math/lp/static_matrix.h"
#include "math/lp/lp_core_solver_base.h"
#include <string>
#include <limits>
#include <set>
#include <algorithm>
#include "util/vector.h"
namespace lp {
template <typename T, typename X>
class lp_dual_core_solver:public lp_core_solver_base<T, X> {
public:
vector<bool> & m_can_enter_basis;
int m_r; // the row of the leaving column
int m_p; // leaving column; that is m_p = m_basis[m_r]
T m_delta; // the offset of the leaving basis variable
int m_sign_of_alpha_r; // see page 27
T m_theta_D;
T m_theta_P;
int m_q;
// todo : replace by a vector later
std::set<unsigned> m_breakpoint_set; // it is F in "Progress in the dual simplex method ..."
std::set<unsigned> m_flipped_boxed;
std::set<unsigned> m_tight_set; // it is the set of all breakpoints that become tight when m_q becomes tight
vector<T> m_a_wave;
vector<T> m_betas; // m_betas[i] is approximately a square of the norm of the i-th row of the reverse of B
T m_harris_tolerance;
std::set<unsigned> m_forbidden_rows;
lp_dual_core_solver(static_matrix<T, X> & A,
vector<bool> & can_enter_basis,
vector<X> & b, // the right side vector
vector<X> & x, // the number of elements in x needs to be at least as large as the number of columns in A
vector<unsigned> & basis,
vector<unsigned> & nbasis,
vector<int> & heading,
vector<T> & costs,
vector<column_type> & column_type_array,
vector<X> & lower_bound_values,
vector<X> & upper_bound_values,
lp_settings & settings,
const column_namer & column_names):
lp_core_solver_base<T, X>(A,
// b,
basis,
nbasis,
heading,
x,
costs,
settings,
column_names,
column_type_array,
lower_bound_values,
upper_bound_values),
m_can_enter_basis(can_enter_basis),
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);
lp_assert(false);
this->init_basic_part_of_basis_heading();
fill_non_basis_with_only_able_to_enter_columns();
}
void init_a_wave_by_zeros();
void fill_non_basis_with_only_able_to_enter_columns() {
auto & nb = this->m_nbasis;
nb.reset();
unsigned j = this->m_n();
while (j--) {
if (this->m_basis_heading[j] >= 0 || !m_can_enter_basis[j]) continue;
nb.push_back(j);
this->m_basis_heading[j] = - static_cast<int>(nb.size());
}
}
void restore_non_basis();
bool update_basis(int entering, int leaving);
void recalculate_xB_and_d();
void recalculate_d();
void init_betas();
void adjust_xb_for_changed_xn_and_init_betas();
void start_with_initial_basis_and_make_it_dual_feasible();
bool done();
T get_edge_steepness_for_lower_bound(unsigned p);
T get_edge_steepness_for_upper_bound(unsigned p);
T pricing_for_row(unsigned i);
void pricing_loop(unsigned number_of_rows_to_try, unsigned offset_in_rows);
bool advance_on_known_p();
int define_sign_of_alpha_r();
bool can_be_breakpoint(unsigned j);
void fill_breakpoint_set();
void DSE_FTran();
T get_delta();
void restore_d();
bool d_is_correct();
void xb_minus_delta_p_pivot_column();
void update_betas();
void apply_flips();
void snap_xN_column_to_bounds(unsigned j);
void snap_xN_to_bounds();
void init_beta_precisely(unsigned i);
void init_betas_precisely();
// step 7 of the algorithm from Progress
bool basis_change_and_update();
void revert_to_previous_basis();
non_basic_column_value_position m_entering_boundary_position;
bool update_basis_and_x_local(int entering, int leaving, X const & tt);
void recover_leaving();
bool problem_is_dual_feasible() const;
bool snap_runaway_nonbasic_column(unsigned);
bool snap_runaway_nonbasic_columns();
unsigned get_number_of_rows_to_try_for_leaving();
void update_a_wave(const T & del, unsigned j) {
this->m_A.add_column_to_vector(del, j, & m_a_wave[0]);
}
bool delta_keeps_the_sign(int initial_delta_sign, const T & delta);
void set_status_to_tentative_dual_unbounded_or_dual_unbounded();
// it is positive if going from low bound to upper bound and negative if going from upper bound to low bound
T signed_span_of_boxed(unsigned j) {
return this->x_is_at_lower_bound(j)? this->bound_span(j): - this->bound_span(j);
}
void add_tight_breakpoints_and_q_to_flipped_set();
T delta_lost_on_flips_of_tight_breakpoints();
bool tight_breakpoinst_are_all_boxed();
T calculate_harris_delta_on_breakpoint_set();
void fill_tight_set_on_harris_delta(const T & harris_delta );
void find_q_on_tight_set();
void find_q_and_tight_set();
void erase_tight_breakpoints_and_q_from_breakpoint_set();
bool ratio_test();
void process_flipped();
void update_d_and_xB();
void calculate_beta_r_precisely();
// see "Progress in the dual simplex method for large scale LP problems: practical dual phase 1 algorithms"
void update_xb_after_bound_flips();
void one_iteration();
void solve();
bool lower_bounds_are_set() const override { return true; }
};
}

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

@ -1,699 +0,0 @@
/*++
Copyright (c) 2017 Microsoft Corporation
Module Name:
<name>
Abstract:
<abstract>
Author:
Lev Nachmanson (levnach)
Revision History:
--*/
#pragma once
#include <algorithm>
#include <string>
#include "util/vector.h"
#include "math/lp/lp_dual_core_solver.h"
namespace lp {
template <typename T, typename X> void lp_dual_core_solver<T, X>::init_a_wave_by_zeros() {
unsigned j = this->m_m();
while (j--) {
m_a_wave[j] = numeric_traits<T>::zero();
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_non_basis() {
auto & nb = this->m_nbasis;
nb.reset();
unsigned j = this->m_n();
while (j--) {
if (this->m_basis_heading[j] >= 0 ) continue;
if (m_can_enter_basis[j]) {
lp_assert(std::find(nb.begin(), nb.end(), j) == nb.end());
nb.push_back(j);
this->m_basis_heading[j] = - static_cast<int>(nb.size());
}
}
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::update_basis(int entering, int leaving) {
// the second argument is the element of the entering column from the pivot row - its value should be equal to the low diagonal element of the bump after all pivoting is done
if (this->m_refactor_counter++ < 200) {
this->m_factorization->replace_column(this->m_ed[this->m_factorization->basis_heading(leaving)], this->m_w);
if (this->m_factorization->get_status() == LU_status::OK) {
this->m_factorization->change_basis(entering, leaving);
return true;
}
}
// need to refactor
this->m_factorization->change_basis(entering, leaving);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings);
this->m_refactor_counter = 0;
if (this->m_factorization->get_status() != LU_status::OK) {
LP_OUT(this->m_settings, "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->total_iterations() << std::endl);
this->m_iters_with_no_cost_growing++;
return false;
}
return true;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::recalculate_xB_and_d() {
this->solve_Ax_eq_b();
recalculate_d();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::recalculate_d() {
lp_assert(false)
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::init_betas() {
// todo : look at page 194 of Progress in the dual simplex algorithm for solving large scale LP problems : techniques for a fast and stable implementation
// the current implementation is not good enough: todo
unsigned i = this->m_m();
while (i--) {
m_betas[i] = 1;
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::adjust_xb_for_changed_xn_and_init_betas() {
this->solve_Ax_eq_b();
init_betas();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::start_with_initial_basis_and_make_it_dual_feasible() {
this->set_non_basic_x_to_correct_bounds(); // It is not an efficient version, see 3.29,
// however this version does not require that m_x is the solution of Ax = 0 beforehand
adjust_xb_for_changed_xn_and_init_betas();
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::done() {
if (this->get_status() == lp_status::OPTIMAL) {
return true;
}
return false; // todo, need to be more cases
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::get_edge_steepness_for_lower_bound(unsigned p) {
lp_assert(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
T del = this->m_x[p] - this->m_lower_bounds[p];
del *= del;
return del / this->m_betas[this->m_basis_heading[p]];
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::get_edge_steepness_for_upper_bound(unsigned p) {
lp_assert(this->m_basis_heading[p] >= 0 && static_cast<unsigned>(this->m_basis_heading[p]) < this->m_m());
T del = this->m_x[p] - this->m_upper_bounds[p];
del *= del;
return del / this->m_betas[this->m_basis_heading[p]];
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::pricing_for_row(unsigned i) {
unsigned p = this->m_basis[i];
switch (this->m_column_types[p]) {
case column_type::fixed:
case column_type::boxed:
if (this->x_below_low_bound(p)) {
T del = get_edge_steepness_for_lower_bound(p);
return del;
}
if (this->x_above_upper_bound(p)) {
T del = get_edge_steepness_for_upper_bound(p);
return del;
}
return numeric_traits<T>::zero();
case column_type::lower_bound:
if (this->x_below_low_bound(p)) {
T del = get_edge_steepness_for_lower_bound(p);
return del;
}
return numeric_traits<T>::zero();
break;
case column_type::upper_bound:
if (this->x_above_upper_bound(p)) {
T del = get_edge_steepness_for_upper_bound(p);
return del;
}
return numeric_traits<T>::zero();
break;
case column_type::free_column:
lp_assert(numeric_traits<T>::is_zero(this->m_d[p]));
return numeric_traits<T>::zero();
default:
lp_unreachable();
}
lp_unreachable();
return numeric_traits<T>::zero();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(unsigned number_of_rows_to_try, unsigned offset_in_rows) {
m_r = -1;
T steepest_edge_max = numeric_traits<T>::zero();
unsigned initial_offset_in_rows = offset_in_rows;
unsigned i = offset_in_rows;
unsigned rows_left = number_of_rows_to_try;
do {
if (m_forbidden_rows.find(i) != m_forbidden_rows.end()) {
if (++i == this->m_m()) {
i = 0;
}
continue;
}
T se = pricing_for_row(i);
if (se > steepest_edge_max) {
steepest_edge_max = se;
m_r = i;
if (rows_left > 0) {
rows_left--;
}
}
if (++i == this->m_m()) {
i = 0;
}
} while (i != initial_offset_in_rows && rows_left);
if (m_r == -1) {
if (this->get_status() != lp_status::UNSTABLE) {
this->set_status(lp_status::OPTIMAL);
}
} else {
m_p = this->m_basis[m_r];
m_delta = get_delta();
if (advance_on_known_p()){
m_forbidden_rows.clear();
return;
}
// failure in advance_on_known_p
if (this->get_status() == lp_status::FLOATING_POINT_ERROR) {
return;
}
this->set_status(lp_status::UNSTABLE);
m_forbidden_rows.insert(m_r);
}
}
// this calculation is needed for the steepest edge update,
// it hijackes m_pivot_row_of_B_1 for this purpose since we will need it anymore to the end of the cycle
template <typename T, typename X> void lp_dual_core_solver<T, X>::DSE_FTran() { // todo, see algorithm 7 from page 35
this->m_factorization->solve_By_for_T_indexed_only(this->m_pivot_row_of_B_1, this->m_settings);
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_known_p() {
return false;
}
template <typename T, typename X> int lp_dual_core_solver<T, X>::define_sign_of_alpha_r() {
switch (this->m_column_types[m_p]) {
case column_type::boxed:
case column_type::fixed:
if (this->x_below_low_bound(m_p)) {
return -1;
}
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lp_unreachable();
case column_type::lower_bound:
if (this->x_below_low_bound(m_p)) {
return -1;
}
lp_unreachable();
case column_type::upper_bound:
if (this->x_above_upper_bound(m_p)) {
return 1;
}
lp_unreachable();
default:
lp_unreachable();
}
lp_unreachable();
return 0;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::can_be_breakpoint(unsigned j) {
if (this->pivot_row_element_is_too_small_for_ratio_test(j)) return false;
switch (this->m_column_types[j]) {
case column_type::lower_bound:
lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_lower_bounds[j]));
return m_sign_of_alpha_r * this->m_pivot_row[j] > 0;
case column_type::upper_bound:
lp_assert(this->m_settings.abs_val_is_smaller_than_harris_tolerance(this->m_x[j] - this->m_upper_bounds[j]));
return m_sign_of_alpha_r * this->m_pivot_row[j] < 0;
case column_type::boxed:
{
bool lower_bound = this->x_is_at_lower_bound(j);
bool grawing = m_sign_of_alpha_r * this->m_pivot_row[j] > 0;
return lower_bound == grawing;
}
case column_type::fixed: // is always dual feasible so we ignore it
return false;
case column_type::free_column:
return true;
default:
return false;
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::fill_breakpoint_set() {
m_breakpoint_set.clear();
for (unsigned j : this->non_basis()) {
if (can_be_breakpoint(j)) {
m_breakpoint_set.insert(j);
}
}
}
// template <typename T, typename X> void lp_dual_core_solver<T, X>::FTran() {
// this->solve_Bd(m_q);
// }
template <typename T, typename X> T lp_dual_core_solver<T, X>::get_delta() {
switch (this->m_column_types[m_p]) {
case column_type::boxed:
if (this->x_below_low_bound(m_p)) {
return this->m_x[m_p] - this->m_lower_bounds[m_p];
}
if (this->x_above_upper_bound(m_p)) {
return this->m_x[m_p] - this->m_upper_bounds[m_p];
}
lp_unreachable();
case column_type::lower_bound:
if (this->x_below_low_bound(m_p)) {
return this->m_x[m_p] - this->m_lower_bounds[m_p];
}
lp_unreachable();
case column_type::upper_bound:
if (this->x_above_upper_bound(m_p)) {
return get_edge_steepness_for_upper_bound(m_p);
}
lp_unreachable();
case column_type::fixed:
return this->m_x[m_p] - this->m_upper_bounds[m_p];
default:
lp_unreachable();
}
lp_unreachable();
return zero_of_type<T>();
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_d() {
this->m_d[m_p] = numeric_traits<T>::zero();
for (auto j : this->non_basis()) {
this->m_d[j] += m_theta_D * this->m_pivot_row[j];
}
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::d_is_correct() {
lp_assert(false);
return true;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::xb_minus_delta_p_pivot_column() {
unsigned i = this->m_m();
while (i--) {
this->m_x[this->m_basis[i]] -= m_theta_P * this->m_ed[i];
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::update_betas() { // page 194 of Progress ... todo - once in a while betas have to be reinitialized
T one_over_arq = numeric_traits<T>::one() / this->m_pivot_row[m_q];
T beta_r = this->m_betas[m_r] = std::max(T(0.0001), (m_betas[m_r] * one_over_arq) * one_over_arq);
T k = -2 * one_over_arq;
unsigned i = this->m_m();
while (i--) {
if (static_cast<int>(i) == m_r) continue;
T a = this->m_ed[i];
m_betas[i] += a * (a * beta_r + k * this->m_pivot_row_of_B_1[i]);
if (m_betas[i] < T(0.0001))
m_betas[i] = T(0.0001);
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::apply_flips() {
for (unsigned j : m_flipped_boxed) {
lp_assert(this->x_is_at_bound(j));
if (this->x_is_at_lower_bound(j)) {
this->m_x[j] = this->m_upper_bounds[j];
} else {
this->m_x[j] = this->m_lower_bounds[j];
}
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::snap_xN_column_to_bounds(unsigned j) {
switch (this->m_column_type[j]) {
case column_type::fixed:
this->m_x[j] = this->m_lower_bounds[j];
break;
case column_type::boxed:
if (this->x_is_at_lower_bound(j)) {
this->m_x[j] = this->m_lower_bounds[j];
} else {
this->m_x[j] = this->m_upper_bounds[j];
}
break;
case column_type::lower_bound:
this->m_x[j] = this->m_lower_bounds[j];
break;
case column_type::upper_bound:
this->m_x[j] = this->m_upper_bounds[j];
break;
case column_type::free_column:
break;
default:
lp_unreachable();
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::snap_xN_to_bounds() {
for (auto j : this->non_basis()) {
snap_xN_column_to_bounds(j);
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::init_beta_precisely(unsigned i) {
vector<T> vec(this->m_m(), numeric_traits<T>::zero());
vec[i] = numeric_traits<T>::one();
this->m_factorization->solve_yB_with_error_check(vec, this->m_basis);
T beta = numeric_traits<T>::zero();
for (T & v : vec) {
beta += v * v;
}
this->m_betas[i] =beta;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::init_betas_precisely() {
unsigned i = this->m_m();
while (i--) {
init_beta_precisely(i);
}
}
// step 7 of the algorithm from Progress
template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_and_update() {
return true;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::recover_leaving() {
switch (m_entering_boundary_position) {
case at_lower_bound:
case at_fixed:
this->m_x[m_q] = this->m_lower_bounds[m_q];
break;
case at_upper_bound:
this->m_x[m_q] = this->m_upper_bounds[m_q];
break;
case free_of_bounds:
this->m_x[m_q] = zero_of_type<X>();
default:
lp_unreachable();
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::revert_to_previous_basis() {
LP_OUT(this->m_settings, "revert to previous basis on ( " << m_p << ", " << m_q << ")" << std::endl);
this->change_basis_unconditionally(m_p, m_q);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_settings);
if (this->m_factorization->get_status() != LU_status::OK) {
this->set_status(lp_status::FLOATING_POINT_ERROR); // complete failure
return;
}
recover_leaving();
if (!this->find_x_by_solving()) {
this->set_status(lp_status::FLOATING_POINT_ERROR);
return;
}
recalculate_xB_and_d();
init_betas_precisely();
}
// returns true if the column has been snapped
template <typename T, typename X> bool lp_dual_core_solver<T, X>::snap_runaway_nonbasic_column(unsigned j) {
switch (this->m_column_types[j]) {
case column_type::fixed:
case column_type::lower_bound:
if (!this->x_is_at_lower_bound(j)) {
this->m_x[j] = this->m_lower_bounds[j];
return true;
}
break;
case column_type::boxed:
{
bool closer_to_lower_bound = abs(this->m_lower_bounds[j] - this->m_x[j]) < abs(this->m_upper_bounds[j] - this->m_x[j]);
if (closer_to_lower_bound) {
if (!this->x_is_at_lower_bound(j)) {
this->m_x[j] = this->m_lower_bounds[j];
return true;
}
} else {
if (!this->x_is_at_upper_bound(j)) {
this->m_x[j] = this->m_lower_bounds[j];
return true;
}
}
}
break;
case column_type::upper_bound:
if (!this->x_is_at_upper_bound(j)) {
this->m_x[j] = this->m_upper_bounds[j];
return true;
}
break;
default:
break;
}
return false;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() const {
for (unsigned j : this->non_basis()){
if (!this->column_is_dual_feasible(j)) {
return false;
}
}
return true;
}
template <typename T, typename X> unsigned lp_dual_core_solver<T, X>::get_number_of_rows_to_try_for_leaving() {
unsigned s = this->m_m();
if (this->m_m() > 300) {
s = (unsigned)((s / 100.0) * this->m_settings.percent_of_entering_to_check);
}
return this->m_settings.random_next() % s + 1;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::delta_keeps_the_sign(int initial_delta_sign, const T & delta) {
if (numeric_traits<T>::precise())
return ((delta > numeric_traits<T>::zero()) && (initial_delta_sign == 1)) ||
((delta < numeric_traits<T>::zero()) && (initial_delta_sign == -1));
double del = numeric_traits<T>::get_double(delta);
return ( (del > this->m_settings.zero_tolerance) && (initial_delta_sign == 1)) ||
((del < - this->m_settings.zero_tolerance) && (initial_delta_sign == -1));
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::set_status_to_tentative_dual_unbounded_or_dual_unbounded() {
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
this->set_status(lp_status::DUAL_UNBOUNDED);
} else {
this->set_status(lp_status::TENTATIVE_DUAL_UNBOUNDED);
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::add_tight_breakpoints_and_q_to_flipped_set() {
m_flipped_boxed.insert(m_q);
for (auto j : m_tight_set) {
m_flipped_boxed.insert(j);
}
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::delta_lost_on_flips_of_tight_breakpoints() {
T ret = abs(this->bound_span(m_q) * this->m_pivot_row[m_q]);
for (auto j : m_tight_set) {
ret += abs(this->bound_span(j) * this->m_pivot_row[j]);
}
return ret;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::tight_breakpoinst_are_all_boxed() {
if (this->m_column_types[m_q] != column_type::boxed) return false;
for (auto j : m_tight_set) {
if (this->m_column_types[j] != column_type::boxed) return false;
}
return true;
}
template <typename T, typename X> T lp_dual_core_solver<T, X>::calculate_harris_delta_on_breakpoint_set() {
bool first_time = true;
T ret = zero_of_type<T>();
lp_assert(m_breakpoint_set.size() > 0);
for (auto j : m_breakpoint_set) {
T t;
if (this->x_is_at_lower_bound(j)) {
t = abs((std::max(this->m_d[j], numeric_traits<T>::zero()) + m_harris_tolerance) / this->m_pivot_row[j]);
} else {
t = abs((std::min(this->m_d[j], numeric_traits<T>::zero()) - m_harris_tolerance) / this->m_pivot_row[j]);
}
if (first_time) {
ret = t;
first_time = false;
} else if (t < ret) {
ret = t;
}
}
return ret;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::fill_tight_set_on_harris_delta(const T & harris_delta ){
m_tight_set.clear();
for (auto j : m_breakpoint_set) {
if (this->x_is_at_lower_bound(j)) {
if (abs(std::max(this->m_d[j], numeric_traits<T>::zero()) / this->m_pivot_row[j]) <= harris_delta){
m_tight_set.insert(j);
}
} else {
if (abs(std::min(this->m_d[j], numeric_traits<T>::zero() ) / this->m_pivot_row[j]) <= harris_delta){
m_tight_set.insert(j);
}
}
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::find_q_on_tight_set() {
m_q = -1;
T max_pivot;
for (auto j : m_tight_set) {
T r = abs(this->m_pivot_row[j]);
if (m_q != -1) {
if (r > max_pivot) {
max_pivot = r;
m_q = j;
}
} else {
max_pivot = r;
m_q = j;
}
}
m_tight_set.erase(m_q);
lp_assert(m_q != -1);
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::find_q_and_tight_set() {
T harris_del = calculate_harris_delta_on_breakpoint_set();
fill_tight_set_on_harris_delta(harris_del);
find_q_on_tight_set();
m_entering_boundary_position = this->get_non_basic_column_value_position(m_q);
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::erase_tight_breakpoints_and_q_from_breakpoint_set() {
m_breakpoint_set.erase(m_q);
for (auto j : m_tight_set) {
m_breakpoint_set.erase(j);
}
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::ratio_test() {
m_sign_of_alpha_r = define_sign_of_alpha_r();
fill_breakpoint_set();
m_flipped_boxed.clear();
int initial_delta_sign = m_delta >= numeric_traits<T>::zero()? 1: -1;
do {
if (m_breakpoint_set.empty()) {
set_status_to_tentative_dual_unbounded_or_dual_unbounded();
return false;
}
this->set_status(lp_status::FEASIBLE);
find_q_and_tight_set();
if (!tight_breakpoinst_are_all_boxed()) break;
T del = m_delta - delta_lost_on_flips_of_tight_breakpoints() * initial_delta_sign;
if (!delta_keeps_the_sign(initial_delta_sign, del)) break;
if (m_tight_set.size() + 1 == m_breakpoint_set.size()) {
break; // deciding not to flip since we might get stuck without finding m_q, the column entering the basis
}
// we can flip m_q together with the tight set and look for another breakpoint candidate for m_q and another tight set
add_tight_breakpoints_and_q_to_flipped_set();
m_delta = del;
erase_tight_breakpoints_and_q_from_breakpoint_set();
} while (true);
m_theta_D = this->m_d[m_q] / this->m_pivot_row[m_q];
return true;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::process_flipped() {
init_a_wave_by_zeros();
for (auto j : m_flipped_boxed) {
update_a_wave(signed_span_of_boxed(j), j);
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::update_d_and_xB() {
for (auto j : this->non_basis()) {
this->m_d[j] -= m_theta_D * this->m_pivot_row[j];
}
this->m_d[m_p] = - m_theta_D;
if (!m_flipped_boxed.empty()) {
process_flipped();
update_xb_after_bound_flips();
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::calculate_beta_r_precisely() {
T t = numeric_traits<T>::zero();
unsigned i = this->m_m();
while (i--) {
T b = this->m_pivot_row_of_B_1[i];
t += b * b;
}
m_betas[m_r] = t;
}
// see "Progress in the dual simplex method for large scale LP problems: practical dual phase 1 algorithms"
template <typename T, typename X> void lp_dual_core_solver<T, X>::update_xb_after_bound_flips() {
this->m_factorization->solve_By(m_a_wave);
unsigned i = this->m_m();
while (i--) {
this->m_x[this->m_basis[i]] -= m_a_wave[i];
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::one_iteration() {
unsigned number_of_rows_to_try = get_number_of_rows_to_try_for_leaving();
unsigned offset_in_rows = this->m_settings.random_next() % this->m_m();
if (this->get_status() == lp_status::TENTATIVE_DUAL_UNBOUNDED) {
number_of_rows_to_try = this->m_m();
} else {
this->set_status(lp_status::FEASIBLE);
}
pricing_loop(number_of_rows_to_try, offset_in_rows);
lp_assert(problem_is_dual_feasible());
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::solve() { // see the page 35
lp_assert(d_is_correct());
lp_assert(problem_is_dual_feasible());
lp_assert(this->basis_heading_is_correct());
//this->set_total_iterations(0);
this->iters_with_no_cost_growing() = 0;
do {
if (this->print_statistics_with_iterations_and_nonzeroes_and_cost_and_check_that_the_time_is_over("", *this->m_settings.get_message_ostream())){
return;
}
one_iteration();
} while (this->get_status() != lp_status::FLOATING_POINT_ERROR && this->get_status() != lp_status::DUAL_UNBOUNDED && this->get_status() != lp_status::OPTIMAL &&
this->iters_with_no_cost_growing() <= this->m_settings.max_number_of_iterations_with_no_improvements
);
}
}

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

@ -23,7 +23,6 @@ Revision History:
#include <utility>
#include "math/lp/static_matrix_def.h"
#include "math/lp/lp_core_solver_base.h"
#include "math/lp/lp_dual_core_solver.h"
#include "math/lp/lp_primal_core_solver.h"
#include "math/lp/scaler.h"
#include "math/lp/lar_solver.h"