michael/lean2
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

dev(lp): fix bugs in the dual simplex high loop

Browse source 
Signed-off-by: Lev Nachmanson <levnach@microsoft.com>
This commit is contained in:
Lev Nachmanson 2016-02-01 15:53:01 -08:00 committed by Leonardo de Moura
parent fbe4f56aea
commit 9482d508df
9 changed files with 210 additions and 113 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

41
src/util/lp/lp_core_solver_base.cpp
View file

@ -198,11 +198,11 @@ A_mult_x_is_off() {
X eps = feps * (one + T(0.1) * abs(m_b[i]));
if (delta > eps) {
std::cout << "x is off (";
std::cout << "m_b[" << i << "] = " << m_b[i] << " ";
std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
std::cout << "delta = " << delta << ' ';
std::cout << "iters = " << m_total_iterations << ")" << std::endl;
// std::cout << "x is off (";
// std::cout << "m_b[" << i << "] = " << m_b[i] << " ";
// std::cout << "left side = " << m_A.dot_product_with_row(i, m_x) << ' ';
// std::cout << "delta = " << delta << ' ';
// std::cout << "iters = " << m_total_iterations << ")" << std::endl;
return true;
}
}
@ -333,7 +333,7 @@ set_non_basic_x_to_correct_bounds() {
}
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
column_is_dual_feasible(unsigned j) {
column_is_dual_feasible(unsigned j) const {
switch (m_column_type[j]) {
case fixed:
case boxed:
@ -353,7 +353,7 @@ column_is_dual_feasible(unsigned j) {
}
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
d_is_not_negative(unsigned j) {
d_is_not_negative(unsigned j) const {
if (numeric_traits<T>::precise()) {
return m_d[j] >= numeric_traits<T>::zero();
}
@ -361,7 +361,7 @@ d_is_not_negative(unsigned j) {
}
template <typename T, typename X> bool lp_core_solver_base<T, X>::
d_is_not_positive(unsigned j) {
d_is_not_positive(unsigned j) const {
if (numeric_traits<T>::precise()) {
return m_d[j] <= numeric_traits<T>::zero();
}
@ -686,4 +686,29 @@ get_non_basic_column_value_position(unsigned j) {
lean_unreachable();
}
}
template <typename T, typename X> void lp_core_solver_base<T, X>::init_lu() {
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
this->m_refactor_counter = 0;
}
template <typename T, typename X> int lp_core_solver_base<T, X>::pivots_in_column_and_row_are_different(int entering, int leaving) const {
const T & column_p = this->m_ed[this->m_basis_heading[leaving]];
const T & row_p = this->m_pivot_row[entering];
if (is_zero(column_p) || is_zero(row_p)) return true; // pivots cannot be zero
// the pivots have to have the same sign
if (column_p < 0) {
if (row_p > 0)
return 2;
} else { // column_p > 0
if (row_p < 0)
return 2;
}
T diff_normalized = abs((column_p - row_p) / (numeric_traits<T>::one() + abs(row_p)));
if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10)))
return 1;
return 0;
}
}

20
src/util/lp/lp_core_solver_base.h
View file

@ -163,6 +163,15 @@ public:
return below_bound(m_x[p], m_low_bound_values[p]);
}
bool x_above_low_bound(unsigned p) {
return above_bound(m_x[p], m_low_bound_values[p]);
}
bool x_below_upper_bound(unsigned p) {
return below_bound(m_x[p], m_upper_bound_values[p]);
}
bool x_above_upper_bound(unsigned p) {
return above_bound(m_x[p], m_upper_bound_values[p]);
}
@ -177,10 +186,11 @@ public:
return x_is_at_low_bound(j) || x_is_at_upper_bound(j);
}
bool column_is_dual_feasible(unsigned j);
bool d_is_not_negative(unsigned j);
bool column_is_dual_feasible(unsigned j) const;
bool d_is_not_positive(unsigned j);
bool d_is_not_negative(unsigned j) const;
bool d_is_not_positive(unsigned j) const;
bool time_is_over();
@ -191,7 +201,6 @@ public:
bool update_basis_and_x(int entering, int leaving, X const & tt);
void init_basis_heading();
bool basis_has_no_doubles();
@ -244,5 +253,8 @@ public:
void init_reduced_costs_for_one_iteration();
non_basic_column_value_position get_non_basic_column_value_position(unsigned j);
void init_lu();
int pivots_in_column_and_row_are_different(int entering, int leaving) const;
};
}

8
src/util/lp/lp_core_solver_base_instances.cpp
View file

@ -9,7 +9,7 @@ template bool lean::lp_core_solver_base<double, double>::A_mult_x_is_off();
template bool lean::lp_core_solver_base<double, double>::basis_heading_is_correct();
template void lean::lp_core_solver_base<double, double>::calculate_pivot_row_of_B_1(unsigned int);
template void lean::lp_core_solver_base<double, double>::calculate_pivot_row_when_pivot_row_of_B1_is_ready();
template bool lean::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int);
template bool lean::lp_core_solver_base<double, double>::column_is_dual_feasible(unsigned int) const;
template void lean::lp_core_solver_base<double, double>::fill_reduced_costs_from_m_y_by_rows();
template bool lean::lp_core_solver_base<double, double>::find_x_by_solving();
template lean::non_basic_column_value_position lean::lp_core_solver_base<double, double>::get_non_basic_column_value_position(unsigned int);
@ -29,7 +29,7 @@ template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::A_mult_x_is_off()
template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::basis_heading_is_correct();
template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::calculate_pivot_row_of_B_1(unsigned int);
template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::calculate_pivot_row_when_pivot_row_of_B1_is_ready();
template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::column_is_dual_feasible(unsigned int);
template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::column_is_dual_feasible(unsigned int) const;
template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::fill_reduced_costs_from_m_y_by_rows();
template bool lean::lp_core_solver_base<lean::mpq, lean::mpq>::find_x_by_solving();
template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::init_reduced_costs_for_one_iteration();
@ -67,3 +67,7 @@ template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq>
template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::restore_state(lean::mpq*, lean::mpq*);
template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::save_state(lean::mpq*, lean::mpq*);
template void lean::lp_core_solver_base<lean::mpq, lean::numeric_pair<lean::mpq> >::solve_yB(std::vector<lean::mpq, std::allocator<lean::mpq> >&);
template void lean::lp_core_solver_base<double, double>::init_lu();
template void lean::lp_core_solver_base<lean::mpq, lean::mpq>::init_lu();
template int lean::lp_core_solver_base<double, double>::pivots_in_column_and_row_are_different(int, int) const;
template int lean::lp_core_solver_base<lean::mpq, lean::mpq>::pivots_in_column_and_row_are_different(int, int) const;

185
src/util/lp/lp_dual_core_solver.cpp
View file

@ -80,29 +80,22 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::restore_non_ba
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)) {
if (this->m_refactor_counter++ < 200) {
this->m_factorization->replace_column(leaving, this->m_ed[this->m_factorization->basis_heading(leaving)], this->m_w);
if (this->m_factorization->get_status() != LU_status::OK) {
#ifdef LEAN_DEBUG
std::cout << "failed on replace_column( " << leaving << ", " << this->m_ed[this->m_factorization->basis_heading(leaving)] << ") total_iterations = " << this->m_total_iterations << std::endl;
#endif
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
this->m_iters_with_no_cost_growing++;
this->m_status = UNSTABLE;
return false;
}
this->m_factorization->change_basis(entering, leaving);
} else { // 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_non_basic_columns);
if (this->m_factorization->get_status() != LU_status::OK) {
#ifdef LEAN_DEBUG
std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << std::endl;
#endif
this->m_iters_with_no_cost_growing++;
return false;
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_non_basic_columns);
this->m_refactor_counter = 0;
if (this->m_factorization->get_status() != LU_status::OK) {
std::cout << "failing refactor for entering = " << entering << ", leaving = " << leaving << " total_iterations = " << this->m_total_iterations << std::endl;
this->m_iters_with_no_cost_growing++;
return false;
}
return true;
}
@ -204,19 +197,24 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
unsigned i = offset_in_rows;
unsigned rows_left = number_of_rows_to_try;
do {
loop_start:
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;
}
if (rows_left > 0) {
rows_left--;
}
} while (rows_left || (steepest_edge_max < T(1e-5) && i != initial_offset_in_rows));
} while (i != initial_offset_in_rows && (rows_left || steepest_edge_max < T(1e-5)));
if (m_r == -1) {
if (this->m_status != UNSTABLE) {
this->m_status = OPTIMAL;
@ -225,22 +223,24 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::pricing_loop(u
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->m_status == FLOATING_POINT_ERROR) {
return;
}
this->m_status = UNSTABLE;
if (i == initial_offset_in_rows) {
return; // made a full loop
}
m_p = m_r = -1;
steepest_edge_max = numeric_traits<T>::zero();
rows_left = number_of_rows_to_try;
goto loop_start;
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(this->m_pivot_row_of_B_1);
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_known_p() {
if (done()) {
return true;
@ -251,7 +251,15 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::advance_on_kno
return true;
}
calculate_beta_r_precisely();
FTran();
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
lean_unreachable(); // not implemented yet
} else {
lean_assert(pivot_compare_result == 1);
this->init_lu();
}
DSE_FTran();
return basis_change_and_update();
}
@ -315,9 +323,9 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::fill_breakpoin
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::FTran() {
this->solve_Bd(m_q);
}
// 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_type[m_p]) {
@ -451,46 +459,115 @@ template <typename T, typename X> void lp_dual_core_solver<T, X>::init_betas_pre
template <typename T, typename X> bool lp_dual_core_solver<T, X>::basis_change_and_update() {
update_betas();
update_d_and_xB();
m_theta_P = m_delta / this->m_ed[m_r];
xb_minus_delta_p_pivot_column();
// m_theta_P = m_delta / this->m_ed[m_r];
m_theta_P = m_delta / this->m_pivot_row[m_q];
// xb_minus_delta_p_pivot_column();
apply_flips();
this->m_x[m_q] += m_theta_P;
if (!update_basis(m_q, m_p) || this->A_mult_x_is_off() || !problem_is_dual_feasible()) {
if (!this->update_basis_and_x(m_q, m_p, m_theta_P)) {
init_betas_precisely();
return false;
}
if (snap_runaway_nonbasic_column(m_p)) {
if(!this->find_x_by_solving()) {
revert_to_previous_basis();
this->m_iters_with_no_cost_growing++;
return false;
}
}
if (!problem_is_dual_feasible() ) {
// todo : shift the costs!!!!
revert_to_previous_basis();
this->m_iters_with_no_cost_growing++;
return false;
}
lean_assert(d_is_correct());
return true;
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::recover_leaving() {
switch (m_entering_boundary_position) {
case at_low_bound:
case at_fixed:
this->m_x[m_q] = this->m_low_bound_values[m_q];
break;
case at_upper_bound:
this->m_x[m_q] = this->m_upper_bound_values[m_q];
break;
case free_of_bounds:
this->m_x[m_q] = zero_of_type<X>();
default:
lean_unreachable();
}
}
template <typename T, typename X> void lp_dual_core_solver<T, X>::revert_to_previous_basis() {
// std::cout << "recovering basis p = " << m_p << " q = " << m_q << std::endl;
this->m_factorization->change_basis(m_p, m_q);
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
if (this->m_factorization->get_status() != LU_status::OK) {
this->m_status = FLOATING_POINT_ERROR; // complete failure
return;
}
recover_leaving();
if (!this->find_x_by_solving()) {
this->m_status = FLOATING_POINT_ERROR;
return;
}
this->m_x[m_q] -= m_theta_P;
snap_xN_to_bounds();
recalculate_xB_and_d();
if (this->A_mult_x_is_off()) {
this->m_status = FLOATING_POINT_ERROR;
return;
}
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_type[j]) {
case fixed:
case low_bound:
if (! this->x_is_at_low_bound(j)) {
this->m_x[j] = this->m_low_bound_values[j];
return true;
}
break;
case boxed:
{
bool closer_to_low_bound = abs(this->m_low_bound_values[j] - this->m_x[j]) < abs(this->m_upper_bound_values[j] - this->m_x[j]);
if (closer_to_low_bound) {
if (! this->x_is_at_low_bound(j)) {
this->m_x[j] = this->m_low_bound_values[j];
return true;
}
} else {
if (! this->x_is_at_upper_bound(j)) {
this->m_x[j] = this->m_low_bound_values[j];
return true;
}
}
}
break;
case upper_bound:
if (! this->x_is_at_upper_bound(j)) {
this->m_x[j] = this->m_upper_bound_values[j];
return true;
}
break;
default:
break;
}
return false;
}
template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() {
template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dual_feasible() const {
for (unsigned j : non_basis()){
if (!this->column_is_dual_feasible(j)) {
std::cout << "column " << j << " is not dual feasible" << std::endl;
std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
std::cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << std::endl;
std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
// std::cout << "column " << j << " is not dual feasible" << std::endl;
// std::cout << "m_d[" << j << "] = " << this->m_d[j] << std::endl;
// std::cout << "x[" << j << "] = " << this->m_x[j] << std::endl;
// std::cout << "type = " << column_type_to_string(this->m_column_type[j]) << std::endl;
// std::cout << "bounds = " << this->m_low_bound_values[j] << "," << this->m_upper_bound_values[j] << std::endl;
// std::cout << "m_total_iterations = " << this->m_total_iterations << std::endl;
return false;
}
}
@ -500,7 +577,7 @@ template <typename T, typename X> bool lp_dual_core_solver<T, X>::problem_is_dua
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 / this->m_settings.percent_of_entering_to_check * 100);
s = (unsigned)((s / 100.0) * this->m_settings.percent_of_entering_to_check);
}
return my_random() % s + 1;
}

23
src/util/lp/lp_dual_core_solver.h
View file

@ -32,7 +32,7 @@ public:
std::vector<T> m_a_wave;
std::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,
std::vector<bool> & can_enter_basis,
std::vector<X> & b, // the right side std::vector
@ -59,7 +59,7 @@ public:
void recalculate_d();
std::vector<unsigned> & non_basis() { return this->m_factorization->m_non_basic_columns; }
const std::vector<unsigned> & non_basis() const { return this->m_factorization->m_non_basic_columns; }
void init_betas();
@ -85,14 +85,7 @@ public:
void fill_breakpoint_set();
void FTran();
// 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
void DSE_FTran() { // todo, see algorithm 7 from page 35
this->m_factorization->solve_By(this->m_pivot_row_of_B_1);
}
void DSE_FTran();
T get_delta();
void restore_d();
@ -118,8 +111,16 @@ public:
void revert_to_previous_basis();
bool problem_is_dual_feasible();
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) {

2
src/util/lp/lp_dual_core_solver_instances.cpp
View file

@ -13,3 +13,5 @@ template void lean::lp_dual_core_solver<double, double>::start_with_initial_basi
template void lean::lp_dual_core_solver<double, double>::solve();
template void lean::lp_dual_core_solver<lean::mpq, lean::mpq>::restore_non_basis();
template void lean::lp_dual_core_solver<double, double>::restore_non_basis();
template void lean::lp_dual_core_solver<double, double>::revert_to_previous_basis();
template void lean::lp_dual_core_solver<lean::mpq, lean::mpq>::revert_to_previous_basis();

5
src/util/lp/lp_dual_simplex.cpp
View file

@ -120,6 +120,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2()
lean_unreachable();
}
this->m_second_stage_iterations = m_core_solver->m_total_iterations;
this->m_total_iterations = this->m_first_stage_iterations + this->m_second_stage_iterations;
}
template <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros() {
@ -272,7 +273,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_s
this->m_costs[artificial] = numeric_traits<T>::zero();
artificial++;
} else {
// we can put a slack_var into the basis, and atemplate <typename T, typename X> void lp_dual_simplex<T, X>::adding an artificial variable
// we can put a slack_var into the basis, and avoid adding an artificial variable
this->m_basis[row] = slack_var;
this->m_costs[slack_var] = numeric_traits<T>::zero();
}
@ -290,7 +291,7 @@ template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_s
this->m_costs[artificial] = numeric_traits<T>::zero();
artificial++;
} else {
// we can put slack_var into the basis, and atemplate <typename T, typename X> void lp_dual_simplex<T, X>::adding an artificial variable
// we can put slack_var into the basis, and avoid adding an artificial variable
this->m_basis[row] = slack_var;
this->m_costs[slack_var] = numeric_traits<T>::zero();
}

36
src/util/lp/lp_primal_core_solver.cpp
View file

@ -261,11 +261,6 @@ lp_primal_core_solver(static_matrix<T, X> & A,
#endif
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::init_lu() {
init_factorization(this->m_factorization, this->m_A, this->m_basis, this->m_basis_heading, this->m_settings, this->m_non_basic_columns);
this->m_refactor_counter = 0;
}
template <typename T, typename X> bool lp_primal_core_solver<T, X>::initial_x_is_correct() {
std::set<unsigned> basis_set;
for (int i = 0; i < this->m_A.row_count(); i++) {
@ -396,23 +391,6 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::init_run(
init_reduced_costs();
}
template <typename T, typename X> int lp_primal_core_solver<T, X>::pivots_in_column_and_row_are_different(int entering, int leaving) const {
const T & column_p = this->m_ed[this->m_basis_heading[leaving]];
const T & row_p = this->m_pivot_row[entering];
if (is_zero(column_p) || is_zero(row_p)) return true; // pivots cannot be zero
// the pivots have to have the same sign
if (column_p < 0) {
if (row_p > 0)
return 2;
} else { // column_p > 0
if (row_p < 0)
return 2;
}
T diff_normalized = abs((column_p - row_p) / (numeric_traits<T>::one() + abs(row_p)));
if ( !this->m_settings.abs_val_is_smaller_than_harris_tolerance(diff_normalized / T(10)))
return 1;
return 0;
}
template <typename T, typename X> void lp_primal_core_solver<T, X>::calc_working_vector_beta_for_column_norms(){
unsigned i = this->m_m;
@ -429,7 +407,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
lean_assert(!this->A_mult_x_is_off() );
this->update_x(entering, t * m_sign_of_entering_delta);
if (this->A_mult_x_is_off() && !this->find_x_by_solving()) {
init_lu();
this->init_lu();
if (!this->find_x_by_solving()) {
this->restore_x(entering, t * m_sign_of_entering_delta);
m_forbidden_enterings.insert(entering);
@ -447,7 +425,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
unsigned pivot_row = this->m_factorization->basis_heading(leaving);
this->calculate_pivot_row_of_B_1(pivot_row);
this->calculate_pivot_row_when_pivot_row_of_B1_is_ready();
int pivot_compare_result = pivots_in_column_and_row_are_different(entering, leaving);
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
m_forbidden_enterings.insert(entering);
@ -456,7 +434,7 @@ template <typename T, typename X>void lp_primal_core_solver<T, X>::advance_on_en
return;
} else {
lean_assert(pivot_compare_result == 1);
init_lu();
this->init_lu();
}
calc_working_vector_beta_for_column_norms();
if (!this->update_basis_and_x(entering, leaving, t * m_sign_of_entering_delta)) {
@ -491,7 +469,7 @@ template <typename T, typename X> void lp_primal_core_solver<T, X>::advance_o
this->solve_Bd(entering);
int refresh_result = refresh_reduced_cost_at_entering_and_check_that_it_is_off(entering);
if (refresh_result) {
init_lu();
this->init_lu();
init_reduced_costs();
if (refresh_result == 2) {
this->m_iters_with_no_cost_growing++;
@ -565,7 +543,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
case OPTIMAL: // double check that we are at optimum
case INFEASIBLE:
m_forbidden_enterings.clear();
init_lu();
this->init_lu();
lean_assert(this->m_factorization->get_status() == LU_status::OK);
set_current_x_is_feasible();
init_reduced_costs();
@ -577,7 +555,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
break;
case TENTATIVE_UNBOUNDED:
m_forbidden_enterings.clear();
init_lu();
this->init_lu();
lean_assert(this->m_factorization->get_status() == LU_status::OK);
init_reduced_costs();
break;
@ -586,7 +564,7 @@ template <typename T, typename X> unsigned lp_primal_core_solver<T, X>::solve()
case UNSTABLE:
// m_forbidden_enterings.clear();
init_lu();
this->init_lu();
init_reduced_costs();
this->m_status = UNKNOWN;
break;

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

@ -127,7 +127,6 @@ public:
std::vector<X> & upper_bound_values,
lp_settings & settings,
std::unordered_map<unsigned, std::string> const & column_names);
void init_lu();
bool initial_x_is_correct();
#ifdef LEAN_DEBUG
void check_Ax_equal_b();
@ -148,8 +147,6 @@ public:
void init_run();
int pivots_in_column_and_row_are_different(int entering, int leaving) const;
void calc_working_vector_beta_for_column_norms();
void advance_on_entering_and_leaving(int entering, int leaving, X & t);