1 /*++ 2 Copyright (c) 2017 Microsoft Corporation 3 Author: 4 5 Lev Nachmanson (levnach) 6 7 --*/ 8 #pragma once 9 #include "util/vector.h" 10 #include <string> 11 #include <utility> 12 #include "math/lp/lp_core_solver_base.h" 13 #include <algorithm> 14 #include "math/lp/indexed_vector.h" 15 #include "math/lp/binary_heap_priority_queue.h" 16 #include "math/lp/breakpoint.h" 17 #include "math/lp/lp_primal_core_solver.h" 18 #include "math/lp/stacked_vector.h" 19 #include "math/lp/lar_solution_signature.h" 20 #include "util/stacked_value.h" 21 namespace lp { 22 23 class lar_core_solver { 24 // m_sign_of_entering is set to 1 if the entering variable needs 25 // to grow and is set to -1 otherwise 26 int m_sign_of_entering_delta; 27 vector<std::pair<mpq, unsigned>> m_infeasible_linear_combination; 28 int m_infeasible_sum_sign; // todo: get rid of this field 29 vector<numeric_pair<mpq>> m_right_sides_dummy; 30 vector<mpq> m_costs_dummy; 31 vector<double> m_d_right_sides_dummy; 32 vector<double> m_d_costs_dummy; 33 public: 34 stacked_value<simplex_strategy_enum> m_stacked_simplex_strategy; 35 stacked_vector<column_type> m_column_types; 36 // r - solver fields, for rational numbers 37 vector<numeric_pair<mpq>> m_r_x; // the solution 38 stacked_vector<numeric_pair<mpq>> m_r_lower_bounds; 39 stacked_vector<numeric_pair<mpq>> m_r_upper_bounds; 40 static_matrix<mpq, numeric_pair<mpq>> m_r_A; 41 stacked_vector<unsigned> m_r_pushed_basis; 42 vector<unsigned> m_r_basis; 43 vector<unsigned> m_r_nbasis; 44 vector<int> m_r_heading; 45 stacked_vector<unsigned> m_r_columns_nz; 46 stacked_vector<unsigned> m_r_rows_nz; 47 48 // d - solver fields, for doubles 49 vector<double> m_d_x; // the solution in doubles 50 vector<double> m_d_lower_bounds; 51 vector<double> m_d_upper_bounds; 52 static_matrix<double, double> m_d_A; 53 stacked_vector<unsigned> m_d_pushed_basis; 54 vector<unsigned> m_d_basis; 55 vector<unsigned> m_d_nbasis; 56 vector<int> m_d_heading; 57 58 59 lp_primal_core_solver<mpq, numeric_pair<mpq>> m_r_solver; // solver in rational numbers 60 61 lp_primal_core_solver<double, double> m_d_solver; // solver in doubles 62 63 lar_core_solver( 64 lp_settings & settings, 65 const column_namer & column_names 66 ); 67 settings()68 lp_settings & settings() { return m_r_solver.m_settings;} settings()69 const lp_settings & settings() const { return m_r_solver.m_settings;} 70 get_infeasible_sum_sign()71 int get_infeasible_sum_sign() const { return m_infeasible_sum_sign; } 72 get_infeasibility_info(int & inf_sign)73 const vector<std::pair<mpq, unsigned>> & get_infeasibility_info(int & inf_sign) const { 74 inf_sign = m_infeasible_sum_sign; 75 return m_infeasible_linear_combination; 76 } 77 78 void fill_not_improvable_zero_sum_from_inf_row(); 79 get_column_type(unsigned j)80 column_type get_column_type(unsigned j) { return m_column_types[j];} 81 82 void calculate_pivot_row(unsigned i); 83 print_pivot_row(std::ostream & out,unsigned row_index)84 void print_pivot_row(std::ostream & out, unsigned row_index) const { 85 for (unsigned j : m_r_solver.m_pivot_row.m_index) { 86 if (numeric_traits<mpq>::is_pos(m_r_solver.m_pivot_row.m_data[j])) 87 out << "+"; 88 out << m_r_solver.m_pivot_row.m_data[j] << m_r_solver.column_name(j) << " "; 89 } 90 91 out << " +" << m_r_solver.column_name(m_r_solver.m_basis[row_index]) << std::endl; 92 93 for (unsigned j : m_r_solver.m_pivot_row.m_index) { 94 m_r_solver.print_column_bound_info(j, out); 95 } 96 m_r_solver.print_column_bound_info(m_r_solver.m_basis[row_index], out); 97 } 98 99 100 void advance_on_sorted_breakpoints(unsigned entering); 101 102 void change_slope_on_breakpoint(unsigned entering, breakpoint<numeric_pair<mpq>> * b, mpq & slope_at_entering); 103 104 bool row_is_infeasible(unsigned row); 105 106 bool row_is_evidence(unsigned row); 107 108 bool find_evidence_row(); 109 110 void prefix_r(); 111 112 void prefix_d(); 113 m_m()114 unsigned m_m() const { return m_r_A.row_count(); } 115 m_n()116 unsigned m_n() const { return m_r_A.column_count(); } 117 is_tiny()118 bool is_tiny() const { return this->m_m() < 10 && this->m_n() < 20; } 119 is_empty()120 bool is_empty() const { return this->m_m() == 0 && this->m_n() == 0; } 121 122 template <typename L> get_sign(const L & v)123 int get_sign(const L & v) { return v > zero_of_type<L>() ? 1 : (v < zero_of_type<L>() ? -1 : 0); } 124 125 void fill_evidence(unsigned row); 126 127 unsigned get_number_of_non_ints() const; 128 129 void solve(); 130 lower_bounds_are_set()131 bool lower_bounds_are_set() const { return true; } 132 get_pivot_row()133 const indexed_vector<mpq> & get_pivot_row() const { 134 return m_r_solver.m_pivot_row; 135 } 136 137 void fill_not_improvable_zero_sum(); 138 pop_basis(unsigned k)139 void pop_basis(unsigned k) { 140 if (!settings().use_tableau()) { 141 m_r_pushed_basis.pop(k); 142 m_r_basis = m_r_pushed_basis(); 143 m_r_solver.init_basis_heading_and_non_basic_columns_vector(); 144 m_d_pushed_basis.pop(k); 145 m_d_basis = m_d_pushed_basis(); 146 m_d_solver.init_basis_heading_and_non_basic_columns_vector(); 147 } else { 148 m_d_basis = m_r_basis; 149 m_d_nbasis = m_r_nbasis; 150 m_d_heading = m_r_heading; 151 } 152 } 153 push()154 void push() { 155 lp_assert(m_r_solver.basis_heading_is_correct()); 156 lp_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct()); 157 lp_assert(m_column_types.size() == m_r_A.column_count()); 158 m_stacked_simplex_strategy = settings().simplex_strategy(); 159 m_stacked_simplex_strategy.push(); 160 m_column_types.push(); 161 // rational 162 if (!settings().use_tableau()) 163 m_r_A.push(); 164 m_r_lower_bounds.push(); 165 m_r_upper_bounds.push(); 166 if (!settings().use_tableau()) { 167 push_vector(m_r_pushed_basis, m_r_basis); 168 push_vector(m_r_columns_nz, m_r_solver.m_columns_nz); 169 push_vector(m_r_rows_nz, m_r_solver.m_rows_nz); 170 } 171 172 m_d_A.push(); 173 if (!settings().use_tableau()) 174 push_vector(m_d_pushed_basis, m_d_basis); 175 } 176 177 template <typename K> push_vector(stacked_vector<K> & pushed_vector,const vector<K> & vector)178 void push_vector(stacked_vector<K> & pushed_vector, const vector<K> & vector) { 179 lp_assert(pushed_vector.size() <= vector.size()); 180 for (unsigned i = 0; i < vector.size();i++) { 181 if (i == pushed_vector.size()) { 182 pushed_vector.push_back(vector[i]); 183 } else { 184 pushed_vector[i] = vector[i]; 185 } 186 } 187 pushed_vector.push(); 188 } 189 pop_markowitz_counts(unsigned k)190 void pop_markowitz_counts(unsigned k) { 191 m_r_columns_nz.pop(k); 192 m_r_rows_nz.pop(k); 193 m_r_solver.m_columns_nz.resize(m_r_columns_nz.size()); 194 m_r_solver.m_rows_nz.resize(m_r_rows_nz.size()); 195 for (unsigned i = 0; i < m_r_columns_nz.size(); i++) 196 m_r_solver.m_columns_nz[i] = m_r_columns_nz[i]; 197 for (unsigned i = 0; i < m_r_rows_nz.size(); i++) 198 m_r_solver.m_rows_nz[i] = m_r_rows_nz[i]; 199 } 200 201 pop(unsigned k)202 void pop(unsigned k) { 203 // rationals 204 if (!settings().use_tableau()) 205 m_r_A.pop(k); 206 m_r_lower_bounds.pop(k); 207 m_r_upper_bounds.pop(k); 208 m_column_types.pop(k); 209 210 delete m_r_solver.m_factorization; 211 m_r_solver.m_factorization = nullptr; 212 m_r_x.resize(m_r_A.column_count()); 213 m_r_solver.m_costs.resize(m_r_A.column_count()); 214 m_r_solver.m_d.resize(m_r_A.column_count()); 215 if(!settings().use_tableau()) 216 pop_markowitz_counts(k); 217 m_d_A.pop(k); 218 // doubles 219 delete m_d_solver.m_factorization; 220 m_d_solver.m_factorization = nullptr; 221 222 m_d_x.resize(m_d_A.column_count()); 223 pop_basis(k); 224 m_stacked_simplex_strategy.pop(k); 225 settings().simplex_strategy() = m_stacked_simplex_strategy; 226 lp_assert(m_r_solver.basis_heading_is_correct()); 227 lp_assert(!need_to_presolve_with_double_solver() || m_d_solver.basis_heading_is_correct()); 228 } 229 need_to_presolve_with_double_solver()230 bool need_to_presolve_with_double_solver() const { 231 return settings().simplex_strategy() == simplex_strategy_enum::lu; 232 } 233 234 template <typename L> is_zero_vector(const vector<L> & b)235 bool is_zero_vector(const vector<L> & b) { 236 for (const L & m: b) 237 if (!is_zero(m)) return false; 238 return true; 239 } 240 241 update_xj_and_get_delta(unsigned j,non_basic_column_value_position pos_type,numeric_pair<mpq> & delta)242 bool update_xj_and_get_delta(unsigned j, non_basic_column_value_position pos_type, numeric_pair<mpq> & delta) { 243 auto & x = m_r_x[j]; 244 switch (pos_type) { 245 case at_lower_bound: 246 if (x == m_r_solver.m_lower_bounds[j]) 247 return false; 248 delta = m_r_solver.m_lower_bounds[j] - x; 249 m_r_solver.m_x[j] = m_r_solver.m_lower_bounds[j]; 250 break; 251 case at_fixed: 252 case at_upper_bound: 253 if (x == m_r_solver.m_upper_bounds[j]) 254 return false; 255 delta = m_r_solver.m_upper_bounds[j] - x; 256 x = m_r_solver.m_upper_bounds[j]; 257 break; 258 case free_of_bounds: { 259 return false; 260 } 261 case not_at_bound: 262 switch (m_column_types[j]) { 263 case column_type::free_column: 264 return false; 265 case column_type::upper_bound: 266 delta = m_r_solver.m_upper_bounds[j] - x; 267 x = m_r_solver.m_upper_bounds[j]; 268 break; 269 case column_type::lower_bound: 270 delta = m_r_solver.m_lower_bounds[j] - x; 271 x = m_r_solver.m_lower_bounds[j]; 272 break; 273 case column_type::boxed: 274 if (x > m_r_solver.m_upper_bounds[j]) { 275 delta = m_r_solver.m_upper_bounds[j] - x; 276 x += m_r_solver.m_upper_bounds[j]; 277 } else { 278 delta = m_r_solver.m_lower_bounds[j] - x; 279 x = m_r_solver.m_lower_bounds[j]; 280 } 281 break; 282 case column_type::fixed: 283 delta = m_r_solver.m_lower_bounds[j] - x; 284 x = m_r_solver.m_lower_bounds[j]; 285 break; 286 287 default: 288 lp_assert(false); 289 } 290 break; 291 default: 292 lp_unreachable(); 293 } 294 m_r_solver.remove_column_from_inf_set(j); 295 return true; 296 } 297 298 299 prepare_solver_x_with_signature_tableau(const lar_solution_signature & signature)300 void prepare_solver_x_with_signature_tableau(const lar_solution_signature & signature) { 301 lp_assert(m_r_solver.inf_set_is_correct()); 302 for (auto &t : signature) { 303 unsigned j = t.first; 304 if (m_r_heading[j] >= 0) 305 continue; 306 auto pos_type = t.second; 307 numeric_pair<mpq> delta; 308 if (!update_xj_and_get_delta(j, pos_type, delta)) 309 continue; 310 for (const auto & cc : m_r_solver.m_A.m_columns[j]){ 311 unsigned i = cc.var(); 312 unsigned jb = m_r_solver.m_basis[i]; 313 m_r_solver.add_delta_to_x_and_track_feasibility(jb, - delta * m_r_solver.m_A.get_val(cc)); 314 } 315 CASSERT("A_off", m_r_solver.A_mult_x_is_off() == false); 316 } 317 lp_assert(m_r_solver.inf_set_is_correct()); 318 } 319 320 321 template <typename L, typename K> prepare_solver_x_with_signature(const lar_solution_signature & signature,lp_primal_core_solver<L,K> & s)322 void prepare_solver_x_with_signature(const lar_solution_signature & signature, lp_primal_core_solver<L,K> & s) { 323 for (auto &t : signature) { 324 unsigned j = t.first; 325 lp_assert(m_r_heading[j] < 0); 326 auto pos_type = t.second; 327 switch (pos_type) { 328 case at_lower_bound: 329 s.m_x[j] = s.m_lower_bounds[j]; 330 break; 331 case at_fixed: 332 case at_upper_bound: 333 s.m_x[j] = s.m_upper_bounds[j]; 334 break; 335 case free_of_bounds: { 336 s.m_x[j] = zero_of_type<K>(); 337 continue; 338 } 339 case not_at_bound: 340 switch (m_column_types[j]) { 341 case column_type::free_column: 342 lp_assert(false); // unreachable 343 case column_type::upper_bound: 344 s.m_x[j] = s.m_upper_bounds[j]; 345 break; 346 case column_type::lower_bound: 347 s.m_x[j] = s.m_lower_bounds[j]; 348 break; 349 case column_type::boxed: 350 if (settings().random_next() % 2) { 351 s.m_x[j] = s.m_lower_bounds[j]; 352 } else { 353 s.m_x[j] = s.m_upper_bounds[j]; 354 } 355 break; 356 case column_type::fixed: 357 s.m_x[j] = s.m_lower_bounds[j]; 358 break; 359 default: 360 lp_assert(false); 361 } 362 break; 363 default: 364 lp_unreachable(); 365 } 366 } 367 368 lp_assert(is_zero_vector(s.m_b)); 369 s.solve_Ax_eq_b(); 370 } 371 372 template <typename L, typename K> catch_up_in_lu_in_reverse(const vector<unsigned> & trace_of_basis_change,lp_primal_core_solver<L,K> & cs)373 void catch_up_in_lu_in_reverse(const vector<unsigned> & trace_of_basis_change, lp_primal_core_solver<L,K> & cs) { 374 // recover the previous working basis 375 for (unsigned i = trace_of_basis_change.size(); i > 0; i-= 2) { 376 unsigned entering = trace_of_basis_change[i-1]; 377 unsigned leaving = trace_of_basis_change[i-2]; 378 cs.change_basis_unconditionally(entering, leaving); 379 } 380 cs.init_lu(); 381 } 382 383 //basis_heading is the basis heading of the solver owning trace_of_basis_change 384 // here we compact the trace as we go to avoid unnecessary column changes 385 template <typename L, typename K> catch_up_in_lu(const vector<unsigned> & trace_of_basis_change,const vector<int> & basis_heading,lp_primal_core_solver<L,K> & cs)386 void catch_up_in_lu(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading, lp_primal_core_solver<L,K> & cs) { 387 if (cs.m_factorization == nullptr || cs.m_factorization->m_refactor_counter + trace_of_basis_change.size()/2 >= 200) { 388 for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { 389 unsigned entering = trace_of_basis_change[i]; 390 unsigned leaving = trace_of_basis_change[i+1]; 391 cs.change_basis_unconditionally(entering, leaving); 392 } 393 if (cs.m_factorization != nullptr) { 394 delete cs.m_factorization; 395 cs.m_factorization = nullptr; 396 } 397 } else { 398 indexed_vector<L> w(cs.m_A.row_count()); 399 // the queues of delayed indices 400 std::queue<unsigned> entr_q, leav_q; 401 auto * l = cs.m_factorization; 402 lp_assert(l->get_status() == LU_status::OK); 403 for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { 404 unsigned entering = trace_of_basis_change[i]; 405 unsigned leaving = trace_of_basis_change[i+1]; 406 bool good_e = basis_heading[entering] >= 0 && cs.m_basis_heading[entering] < 0; 407 bool good_l = basis_heading[leaving] < 0 && cs.m_basis_heading[leaving] >= 0; 408 if (!good_e && !good_l) continue; 409 if (good_e && !good_l) { 410 while (!leav_q.empty() && cs.m_basis_heading[leav_q.front()] < 0) 411 leav_q.pop(); 412 if (!leav_q.empty()) { 413 leaving = leav_q.front(); 414 leav_q.pop(); 415 } else { 416 entr_q.push(entering); 417 continue; 418 } 419 } else if (!good_e && good_l) { 420 while (!entr_q.empty() && cs.m_basis_heading[entr_q.front()] >= 0) 421 entr_q.pop(); 422 if (!entr_q.empty()) { 423 entering = entr_q.front(); 424 entr_q.pop(); 425 } else { 426 leav_q.push(leaving); 427 continue; 428 } 429 } 430 lp_assert(cs.m_basis_heading[entering] < 0); 431 lp_assert(cs.m_basis_heading[leaving] >= 0); 432 if (l->get_status() == LU_status::OK) { 433 l->prepare_entering(entering, w); // to init vector w 434 l->replace_column(zero_of_type<L>(), w, cs.m_basis_heading[leaving]); 435 } 436 cs.change_basis_unconditionally(entering, leaving); 437 } 438 if (l->get_status() != LU_status::OK) { 439 delete l; 440 cs.m_factorization = nullptr; 441 } 442 } 443 if (cs.m_factorization == nullptr) { 444 if (numeric_traits<L>::precise()) 445 init_factorization(cs.m_factorization, cs.m_A, cs.m_basis, settings()); 446 } 447 } 448 no_r_lu()449 bool no_r_lu() const { 450 return m_r_solver.m_factorization == nullptr || m_r_solver.m_factorization->get_status() == LU_status::Degenerated; 451 } 452 solve_on_signature_tableau(const lar_solution_signature & signature,const vector<unsigned> & changes_of_basis)453 void solve_on_signature_tableau(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) { 454 r_basis_is_OK(); 455 lp_assert(settings().use_tableau()); 456 bool r = catch_up_in_lu_tableau(changes_of_basis, m_d_solver.m_basis_heading); 457 458 if (!r) { // it is the case where m_d_solver gives a degenerated basis 459 prepare_solver_x_with_signature_tableau(signature); // still are going to use the signature partially 460 m_r_solver.find_feasible_solution(); 461 m_d_basis = m_r_basis; 462 m_d_heading = m_r_heading; 463 m_d_nbasis = m_r_nbasis; 464 delete m_d_solver.m_factorization; 465 m_d_solver.m_factorization = nullptr; 466 } else { 467 prepare_solver_x_with_signature_tableau(signature); 468 m_r_solver.start_tracing_basis_changes(); 469 m_r_solver.find_feasible_solution(); 470 if (settings().get_cancel_flag()) 471 return; 472 m_r_solver.stop_tracing_basis_changes(); 473 // and now catch up in the double solver 474 lp_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2); 475 catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver); 476 } 477 lp_assert(r_basis_is_OK()); 478 } 479 adjust_x_of_column(unsigned j)480 bool adjust_x_of_column(unsigned j) { 481 /* 482 if (m_r_solver.m_basis_heading[j] >= 0) { 483 return false; 484 } 485 486 if (m_r_solver.column_is_feasible(j)) { 487 return false; 488 } 489 490 m_r_solver.snap_column_to_bound_tableau(j); 491 lp_assert(m_r_solver.column_is_feasible(j)); 492 m_r_solver.m_inf_set.erase(j); 493 */ 494 lp_assert(false); 495 return true; 496 } 497 498 catch_up_in_lu_tableau(const vector<unsigned> & trace_of_basis_change,const vector<int> & basis_heading)499 bool catch_up_in_lu_tableau(const vector<unsigned> & trace_of_basis_change, const vector<int> & basis_heading) { 500 lp_assert(r_basis_is_OK()); 501 // the queues of delayed indices 502 std::queue<unsigned> entr_q, leav_q; 503 for (unsigned i = 0; i < trace_of_basis_change.size(); i+= 2) { 504 unsigned entering = trace_of_basis_change[i]; 505 unsigned leaving = trace_of_basis_change[i+1]; 506 bool good_e = basis_heading[entering] >= 0 && m_r_solver.m_basis_heading[entering] < 0; 507 bool good_l = basis_heading[leaving] < 0 && m_r_solver.m_basis_heading[leaving] >= 0; 508 if (!good_e && !good_l) continue; 509 if (good_e && !good_l) { 510 while (!leav_q.empty() && m_r_solver.m_basis_heading[leav_q.front()] < 0) 511 leav_q.pop(); 512 if (!leav_q.empty()) { 513 leaving = leav_q.front(); 514 leav_q.pop(); 515 } else { 516 entr_q.push(entering); 517 continue; 518 } 519 } else if (!good_e && good_l) { 520 while (!entr_q.empty() && m_r_solver.m_basis_heading[entr_q.front()] >= 0) 521 entr_q.pop(); 522 if (!entr_q.empty()) { 523 entering = entr_q.front(); 524 entr_q.pop(); 525 } else { 526 leav_q.push(leaving); 527 continue; 528 } 529 } 530 lp_assert(m_r_solver.m_basis_heading[entering] < 0); 531 lp_assert(m_r_solver.m_basis_heading[leaving] >= 0); 532 m_r_solver.change_basis_unconditionally(entering, leaving); 533 if(!m_r_solver.pivot_column_tableau(entering, m_r_solver.m_basis_heading[entering])) { 534 // unroll the last step 535 m_r_solver.change_basis_unconditionally(leaving, entering); 536 #ifdef Z3DEBUG 537 bool t = 538 #endif 539 m_r_solver.pivot_column_tableau(leaving, m_r_solver.m_basis_heading[leaving]); 540 #ifdef Z3DEBUG 541 lp_assert(t); 542 #endif 543 return false; 544 } 545 } 546 lp_assert(r_basis_is_OK()); 547 return true; 548 } 549 550 r_basis_is_OK()551 bool r_basis_is_OK() const { 552 #ifdef Z3DEBUG 553 if (!m_r_solver.m_settings.use_tableau()) 554 return true; 555 for (unsigned j : m_r_solver.m_basis) { 556 lp_assert(m_r_solver.m_A.m_columns[j].size() == 1); 557 } 558 for (unsigned j =0; j < m_r_solver.m_basis_heading.size(); j++) { 559 if (m_r_solver.m_basis_heading[j] >= 0) continue; 560 if (m_r_solver.m_column_types[j] == column_type::fixed) continue; 561 lp_assert(static_cast<unsigned>(- m_r_solver.m_basis_heading[j] - 1) < m_r_solver.m_column_types.size()); 562 lp_assert( m_r_solver.m_basis_heading[j] <= -1); 563 } 564 #endif 565 return true; 566 } 567 solve_on_signature(const lar_solution_signature & signature,const vector<unsigned> & changes_of_basis)568 void solve_on_signature(const lar_solution_signature & signature, const vector<unsigned> & changes_of_basis) { 569 SASSERT(!settings().use_tableau()); 570 if (m_r_solver.m_factorization == nullptr) { 571 for (unsigned j = 0; j < changes_of_basis.size(); j+=2) { 572 unsigned entering = changes_of_basis[j]; 573 unsigned leaving = changes_of_basis[j + 1]; 574 m_r_solver.change_basis_unconditionally(entering, leaving); 575 } 576 init_factorization(m_r_solver.m_factorization, m_r_A, m_r_basis, settings()); 577 } else { 578 catch_up_in_lu(changes_of_basis, m_d_solver.m_basis_heading, m_r_solver); 579 } 580 581 if (no_r_lu()) { // it is the case where m_d_solver gives a degenerated basis, we need to roll back 582 catch_up_in_lu_in_reverse(changes_of_basis, m_r_solver); 583 m_r_solver.find_feasible_solution(); 584 m_d_basis = m_r_basis; 585 m_d_heading = m_r_heading; 586 m_d_nbasis = m_r_nbasis; 587 delete m_d_solver.m_factorization; 588 m_d_solver.m_factorization = nullptr; 589 } else { 590 prepare_solver_x_with_signature(signature, m_r_solver); 591 m_r_solver.start_tracing_basis_changes(); 592 m_r_solver.find_feasible_solution(); 593 if (settings().get_cancel_flag()) 594 return; 595 m_r_solver.stop_tracing_basis_changes(); 596 // and now catch up in the double solver 597 lp_assert(m_r_solver.total_iterations() >= m_r_solver.m_trace_of_basis_change_vector.size() /2); 598 catch_up_in_lu(m_r_solver.m_trace_of_basis_change_vector, m_r_solver.m_basis_heading, m_d_solver); 599 } 600 } 601 create_double_matrix(static_matrix<double,double> & A)602 void create_double_matrix(static_matrix<double, double> & A) { 603 for (unsigned i = 0; i < m_r_A.row_count(); i++) { 604 auto & row = m_r_A.m_rows[i]; 605 for (row_cell<mpq> & c : row) { 606 A.add_new_element(i, c.var(), c.coeff().get_double()); 607 } 608 } 609 } 610 fill_basis_d(vector<unsigned> & basis_d,vector<int> & heading_d,vector<unsigned> & nbasis_d)611 void fill_basis_d( 612 vector<unsigned>& basis_d, 613 vector<int>& heading_d, 614 vector<unsigned>& nbasis_d){ 615 basis_d = m_r_basis; 616 heading_d = m_r_heading; 617 nbasis_d = m_r_nbasis; 618 } 619 620 template <typename L, typename K> extract_signature_from_lp_core_solver(const lp_primal_core_solver<L,K> & solver,lar_solution_signature & signature)621 void extract_signature_from_lp_core_solver(const lp_primal_core_solver<L, K> & solver, lar_solution_signature & signature) { 622 signature.clear(); 623 lp_assert(signature.size() == 0); 624 for (unsigned j = 0; j < solver.m_basis_heading.size(); j++) { 625 if (solver.m_basis_heading[j] < 0) { 626 signature[j] = solver.get_non_basic_column_value_position(j); 627 } 628 } 629 } 630 get_bounds_for_double_solver()631 void get_bounds_for_double_solver() { 632 unsigned n = m_n(); 633 m_d_lower_bounds.resize(n); 634 m_d_upper_bounds.resize(n); 635 double delta = find_delta_for_strict_boxed_bounds().get_double(); 636 if (delta > 0.000001) 637 delta = 0.000001; 638 for (unsigned j = 0; j < n; j++) { 639 if (lower_bound_is_set(j)) { 640 const auto & lb = m_r_solver.m_lower_bounds[j]; 641 m_d_lower_bounds[j] = lb.x.get_double() + delta * lb.y.get_double(); 642 } 643 if (upper_bound_is_set(j)) { 644 const auto & ub = m_r_solver.m_upper_bounds[j]; 645 m_d_upper_bounds[j] = ub.x.get_double() + delta * ub.y.get_double(); 646 lp_assert(!lower_bound_is_set(j) || (m_d_upper_bounds[j] >= m_d_lower_bounds[j])); 647 } 648 } 649 } 650 scale_problem_for_doubles(static_matrix<double,double> & A,vector<double> & lower_bounds,vector<double> & upper_bounds)651 void scale_problem_for_doubles( 652 static_matrix<double, double>& A, 653 vector<double> & lower_bounds, 654 vector<double> & upper_bounds) { 655 vector<double> column_scale_vector; 656 vector<double> right_side_vector(A.column_count()); 657 settings().reps_in_scaler = 5; 658 scaler<double, double > scaler(right_side_vector, 659 A, 660 settings().scaling_minimum, 661 settings().scaling_maximum, 662 column_scale_vector, 663 settings()); 664 if (! scaler.scale()) { 665 // the scale did not succeed, unscaling 666 A.clear(); 667 create_double_matrix(A); 668 } else { 669 for (unsigned j = 0; j < A.column_count(); j++) { 670 if (m_r_solver.column_has_upper_bound(j)) { 671 upper_bounds[j] /= column_scale_vector[j]; 672 } 673 if (m_r_solver.column_has_lower_bound(j)) { 674 lower_bounds[j] /= column_scale_vector[j]; 675 } 676 } 677 } 678 679 } 680 // returns the trace of basis changes find_solution_signature_with_doubles(lar_solution_signature & signature)681 vector<unsigned> find_solution_signature_with_doubles(lar_solution_signature & signature) { 682 if (m_d_solver.m_factorization == nullptr || m_d_solver.m_factorization->get_status() != LU_status::OK) { 683 vector<unsigned> ret; 684 return ret; 685 } 686 get_bounds_for_double_solver(); 687 688 extract_signature_from_lp_core_solver(m_r_solver, signature); 689 prepare_solver_x_with_signature(signature, m_d_solver); 690 m_d_solver.start_tracing_basis_changes(); 691 m_d_solver.find_feasible_solution(); 692 if (settings().get_cancel_flag()) 693 return vector<unsigned>(); 694 695 m_d_solver.stop_tracing_basis_changes(); 696 extract_signature_from_lp_core_solver(m_d_solver, signature); 697 return m_d_solver.m_trace_of_basis_change_vector; 698 } 699 700 lower_bound_is_set(unsigned j)701 bool lower_bound_is_set(unsigned j) const { 702 switch (m_column_types[j]) { 703 case column_type::free_column: 704 case column_type::upper_bound: 705 return false; 706 case column_type::lower_bound: 707 case column_type::boxed: 708 case column_type::fixed: 709 return true; 710 default: 711 lp_assert(false); 712 } 713 return false; 714 } 715 upper_bound_is_set(unsigned j)716 bool upper_bound_is_set(unsigned j) const { 717 switch (m_column_types[j]) { 718 case column_type::free_column: 719 case column_type::lower_bound: 720 return false; 721 case column_type::upper_bound: 722 case column_type::boxed: 723 case column_type::fixed: 724 return true; 725 default: 726 lp_assert(false); 727 } 728 return false; 729 } 730 update_delta(mpq & delta,numeric_pair<mpq> const & l,numeric_pair<mpq> const & u)731 void update_delta(mpq& delta, numeric_pair<mpq> const& l, numeric_pair<mpq> const& u) const { 732 lp_assert(l <= u); 733 if (l.x < u.x && l.y > u.y) { 734 mpq delta1 = (u.x - l.x) / (l.y - u.y); 735 if (delta1 < delta) { 736 delta = delta1; 737 } 738 } 739 lp_assert(l.x + delta * l.y <= u.x + delta * u.y); 740 } 741 742 find_delta_for_strict_boxed_bounds()743 mpq find_delta_for_strict_boxed_bounds() const{ 744 mpq delta = numeric_traits<mpq>::one(); 745 for (unsigned j = 0; j < m_r_A.column_count(); j++ ) { 746 if (m_column_types()[j] != column_type::boxed) 747 continue; 748 update_delta(delta, m_r_lower_bounds[j], m_r_upper_bounds[j]); 749 } 750 return delta; 751 } 752 753 find_delta_for_strict_bounds(const mpq & initial_delta)754 mpq find_delta_for_strict_bounds(const mpq & initial_delta) const{ 755 mpq delta = initial_delta; 756 for (unsigned j = 0; j < m_r_A.column_count(); j++ ) { 757 if (lower_bound_is_set(j)) 758 update_delta(delta, m_r_lower_bounds[j], m_r_x[j]); 759 if (upper_bound_is_set(j)) 760 update_delta(delta, m_r_x[j], m_r_upper_bounds[j]); 761 } 762 return delta; 763 } 764 init_column_row_nz_for_r_solver()765 void init_column_row_nz_for_r_solver() { 766 m_r_solver.init_column_row_non_zeroes(); 767 } 768 column_is_fixed(unsigned j)769 bool column_is_fixed(unsigned j) const { 770 return m_column_types()[j] == column_type::fixed || 771 ( m_column_types()[j] == column_type::boxed && 772 m_r_solver.m_lower_bounds[j] == m_r_solver.m_upper_bounds[j]); 773 } 774 column_is_free(unsigned j)775 bool column_is_free(unsigned j) const { 776 return m_column_types()[j] == column_type::free_column; 777 } 778 779 lower_bound(unsigned j)780 const impq & lower_bound(unsigned j) const { 781 lp_assert(m_column_types()[j] == column_type::fixed || 782 m_column_types()[j] == column_type::boxed || 783 m_column_types()[j] == column_type::lower_bound); 784 return m_r_lower_bounds[j]; 785 } 786 upper_bound(unsigned j)787 const impq & upper_bound(unsigned j) const { 788 lp_assert(m_column_types()[j] == column_type::fixed || 789 m_column_types()[j] == column_type::boxed || 790 m_column_types()[j] == column_type::upper_bound); 791 return m_r_upper_bounds[j]; 792 } 793 794 column_is_bounded(unsigned j)795 bool column_is_bounded(unsigned j) const { 796 switch(m_column_types()[j]) { 797 case column_type::fixed: 798 case column_type::boxed: 799 return true; 800 default: 801 return false; 802 } 803 } 804 r_basis()805 const vector<unsigned>& r_basis() const { return m_r_basis; } r_nbasis()806 const vector<unsigned>& r_nbasis() const { return m_r_nbasis; } 807 }; 808 } 809