1 /*++
2 Copyright (c) 2017 Microsoft Corporation
3
4 Module Name:
5
6 <name>
7
8 Abstract:
9
10 <abstract>
11
12 Author:
13
14 Lev Nachmanson (levnach)
15
16 Revision History:
17
18
19 --*/
20 #include "math/lp/lp_dual_simplex.h"
21 namespace lp{
22
decide_on_status_after_stage1()23 template <typename T, typename X> void lp_dual_simplex<T, X>::decide_on_status_after_stage1() {
24 switch (m_core_solver->get_status()) {
25 case lp_status::OPTIMAL:
26 if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
27 this->m_status = lp_status::FEASIBLE;
28 } else {
29 this->m_status = lp_status::UNBOUNDED;
30 }
31 break;
32 case lp_status::DUAL_UNBOUNDED:
33 lp_unreachable();
34 case lp_status::ITERATIONS_EXHAUSTED:
35 this->m_status = lp_status::ITERATIONS_EXHAUSTED;
36 break;
37 case lp_status::TIME_EXHAUSTED:
38 this->m_status = lp_status::TIME_EXHAUSTED;
39 break;
40 case lp_status::FLOATING_POINT_ERROR:
41 this->m_status = lp_status::FLOATING_POINT_ERROR;
42 break;
43 default:
44 lp_unreachable();
45 }
46 }
47
fix_logical_for_stage2(unsigned j)48 template <typename T, typename X> void lp_dual_simplex<T, X>::fix_logical_for_stage2(unsigned j) {
49 lp_assert(j >= this->number_of_core_structurals());
50 switch (m_column_types_of_logicals[j - this->number_of_core_structurals()]) {
51 case column_type::lower_bound:
52 m_lower_bounds[j] = numeric_traits<T>::zero();
53 m_column_types_of_core_solver[j] = column_type::lower_bound;
54 m_can_enter_basis[j] = true;
55 break;
56 case column_type::fixed:
57 this->m_upper_bounds[j] = m_lower_bounds[j] = numeric_traits<T>::zero();
58 m_column_types_of_core_solver[j] = column_type::fixed;
59 m_can_enter_basis[j] = false;
60 break;
61 default:
62 lp_unreachable();
63 }
64 }
65
fix_structural_for_stage2(unsigned j)66 template <typename T, typename X> void lp_dual_simplex<T, X>::fix_structural_for_stage2(unsigned j) {
67 column_info<T> * ci = this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]];
68 switch (ci->get_column_type()) {
69 case column_type::lower_bound:
70 m_lower_bounds[j] = numeric_traits<T>::zero();
71 m_column_types_of_core_solver[j] = column_type::lower_bound;
72 m_can_enter_basis[j] = true;
73 break;
74 case column_type::fixed:
75 case column_type::upper_bound:
76 lp_unreachable();
77 case column_type::boxed:
78 this->m_upper_bounds[j] = ci->get_adjusted_upper_bound() / this->m_column_scale[j];
79 m_lower_bounds[j] = numeric_traits<T>::zero();
80 m_column_types_of_core_solver[j] = column_type::boxed;
81 m_can_enter_basis[j] = true;
82 break;
83 case column_type::free_column:
84 m_can_enter_basis[j] = true;
85 m_column_types_of_core_solver[j] = column_type::free_column;
86 break;
87 default:
88 lp_unreachable();
89 }
90 // T cost_was = this->m_costs[j];
91 this->set_scaled_cost(j);
92 }
93
unmark_boxed_and_fixed_columns_and_fix_structural_costs()94 template <typename T, typename X> void lp_dual_simplex<T, X>::unmark_boxed_and_fixed_columns_and_fix_structural_costs() {
95 unsigned j = this->m_A->column_count();
96 while (j-- > this->number_of_core_structurals()) {
97 fix_logical_for_stage2(j);
98 }
99 j = this->number_of_core_structurals();
100 while (j--) {
101 fix_structural_for_stage2(j);
102 }
103 }
104
restore_right_sides()105 template <typename T, typename X> void lp_dual_simplex<T, X>::restore_right_sides() {
106 unsigned i = this->m_A->row_count();
107 while (i--) {
108 this->m_b[i] = m_b_copy[i];
109 }
110 }
111
solve_for_stage2()112 template <typename T, typename X> void lp_dual_simplex<T, X>::solve_for_stage2() {
113 m_core_solver->restore_non_basis();
114 m_core_solver->solve_yB(m_core_solver->m_y);
115 m_core_solver->fill_reduced_costs_from_m_y_by_rows();
116 m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
117 m_core_solver->set_status(lp_status::FEASIBLE);
118 m_core_solver->solve();
119 switch (m_core_solver->get_status()) {
120 case lp_status::OPTIMAL:
121 this->m_status = lp_status::OPTIMAL;
122 break;
123 case lp_status::DUAL_UNBOUNDED:
124 this->m_status = lp_status::INFEASIBLE;
125 break;
126 case lp_status::TIME_EXHAUSTED:
127 this->m_status = lp_status::TIME_EXHAUSTED;
128 break;
129 case lp_status::FLOATING_POINT_ERROR:
130 this->m_status = lp_status::FLOATING_POINT_ERROR;
131 break;
132 default:
133 lp_unreachable();
134 }
135 this->m_second_stage_iterations = m_core_solver->total_iterations();
136 this->m_total_iterations = (this->m_first_stage_iterations + this->m_second_stage_iterations);
137 }
138
fill_x_with_zeros()139 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_x_with_zeros() {
140 unsigned j = this->m_A->column_count();
141 while (j--) {
142 this->m_x[j] = numeric_traits<T>::zero();
143 }
144 }
145
stage1()146 template <typename T, typename X> void lp_dual_simplex<T, X>::stage1() {
147 lp_assert(m_core_solver == nullptr);
148 this->m_x.resize(this->m_A->column_count(), numeric_traits<T>::zero());
149 if (this->m_settings.get_message_ostream() != nullptr)
150 this->print_statistics_on_A(*this->m_settings.get_message_ostream());
151 m_core_solver = new lp_dual_core_solver<T, X>(
152 *this->m_A,
153 m_can_enter_basis,
154 this->m_b, // the right side vector
155 this->m_x,
156 this->m_basis,
157 this->m_nbasis,
158 this->m_heading,
159 this->m_costs,
160 this->m_column_types_of_core_solver,
161 this->m_lower_bounds,
162 this->m_upper_bounds,
163 this->m_settings,
164 *this);
165 m_core_solver->fill_reduced_costs_from_m_y_by_rows();
166 m_core_solver->start_with_initial_basis_and_make_it_dual_feasible();
167 if (this->m_settings.abs_val_is_smaller_than_artificial_tolerance(m_core_solver->get_cost())) {
168 // skipping stage 1
169 m_core_solver->set_status(lp_status::OPTIMAL);
170 m_core_solver->set_total_iterations(0);
171 } else {
172 m_core_solver->solve();
173 }
174 decide_on_status_after_stage1();
175 this->m_first_stage_iterations = m_core_solver->total_iterations();
176 }
177
stage2()178 template <typename T, typename X> void lp_dual_simplex<T, X>::stage2() {
179 unmark_boxed_and_fixed_columns_and_fix_structural_costs();
180 restore_right_sides();
181 solve_for_stage2();
182 }
183
fill_first_stage_solver_fields()184 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields() {
185 unsigned slack_var = this->number_of_core_structurals();
186 unsigned artificial = this->number_of_core_structurals() + this->m_slacks;
187
188 for (unsigned row = 0; row < this->row_count(); row++) {
189 fill_first_stage_solver_fields_for_row_slack_and_artificial(row, slack_var, artificial);
190 }
191 fill_costs_and_bounds_and_column_types_for_the_first_stage_solver();
192 }
193
get_column_type(unsigned j)194 template <typename T, typename X> column_type lp_dual_simplex<T, X>::get_column_type(unsigned j) {
195 lp_assert(j < this->m_A->column_count());
196 if (j >= this->number_of_core_structurals()) {
197 return m_column_types_of_logicals[j - this->number_of_core_structurals()];
198 }
199 return this->m_map_from_var_index_to_column_info[this->m_core_solver_columns_to_external_columns[j]]->get_column_type();
200 }
201
fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j)202 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(unsigned j) {
203 // see 4.7 in the dissertation of Achim Koberstein
204 lp_assert(this->m_core_solver_columns_to_external_columns.find(j) !=
205 this->m_core_solver_columns_to_external_columns.end());
206
207 T free_bound = T(1e4); // see 4.8
208 unsigned jj = this->m_core_solver_columns_to_external_columns[j];
209 lp_assert(this->m_map_from_var_index_to_column_info.find(jj) != this->m_map_from_var_index_to_column_info.end());
210 column_info<T> * ci = this->m_map_from_var_index_to_column_info[jj];
211 switch (ci->get_column_type()) {
212 case column_type::upper_bound: {
213 std::stringstream s;
214 s << "unexpected bound type " << j << " "
215 << column_type_to_string(get_column_type(j));
216 throw_exception(s.str());
217 break;
218 }
219 case column_type::lower_bound: {
220 m_can_enter_basis[j] = true;
221 this->set_scaled_cost(j);
222 this->m_lower_bounds[j] = numeric_traits<T>::zero();
223 this->m_upper_bounds[j] =numeric_traits<T>::one();
224 break;
225 }
226 case column_type::free_column: {
227 m_can_enter_basis[j] = true;
228 this->set_scaled_cost(j);
229 this->m_upper_bounds[j] = free_bound;
230 this->m_lower_bounds[j] = -free_bound;
231 break;
232 }
233 case column_type::boxed:
234 m_can_enter_basis[j] = false;
235 this->m_costs[j] = numeric_traits<T>::zero();
236 this->m_upper_bounds[j] = this->m_lower_bounds[j] = numeric_traits<T>::zero(); // is it needed?
237 break;
238 default:
239 lp_unreachable();
240 }
241 m_column_types_of_core_solver[j] = column_type::boxed;
242 }
243
fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j)244 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(unsigned j) {
245 this->m_costs[j] = 0;
246 lp_assert(get_column_type(j) != column_type::upper_bound);
247 if ((m_can_enter_basis[j] = (get_column_type(j) == column_type::lower_bound))) {
248 m_column_types_of_core_solver[j] = column_type::boxed;
249 this->m_lower_bounds[j] = numeric_traits<T>::zero();
250 this->m_upper_bounds[j] = numeric_traits<T>::one();
251 } else {
252 m_column_types_of_core_solver[j] = column_type::fixed;
253 this->m_lower_bounds[j] = numeric_traits<T>::zero();
254 this->m_upper_bounds[j] = numeric_traits<T>::zero();
255 }
256 }
257
fill_costs_and_bounds_and_column_types_for_the_first_stage_solver()258 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_costs_and_bounds_and_column_types_for_the_first_stage_solver() {
259 unsigned j = this->m_A->column_count();
260 while (j-- > this->number_of_core_structurals()) { // go over logicals here
261 fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_logical_column(j);
262 }
263 j = this->number_of_core_structurals();
264 while (j--) {
265 fill_costs_bounds_types_and_can_enter_basis_for_the_first_stage_solver_structural_column(j);
266 }
267 }
268
fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,unsigned & slack_var,unsigned & artificial)269 template <typename T, typename X> void lp_dual_simplex<T, X>::fill_first_stage_solver_fields_for_row_slack_and_artificial(unsigned row,
270 unsigned & slack_var,
271 unsigned & artificial) {
272 lp_assert(row < this->row_count());
273 auto & constraint = this->m_constraints[this->m_core_solver_rows_to_external_rows[row]];
274 // we need to bring the program to the form Ax = b
275 T rs = this->m_b[row];
276 switch (constraint.m_relation) {
277 case Equal: // no slack variable here
278 set_type_for_logical(artificial, column_type::fixed);
279 this->m_basis[row] = artificial;
280 this->m_costs[artificial] = numeric_traits<T>::zero();
281 (*this->m_A)(row, artificial) = numeric_traits<T>::one();
282 artificial++;
283 break;
284
285 case Greater_or_equal:
286 set_type_for_logical(slack_var, column_type::lower_bound);
287 (*this->m_A)(row, slack_var) = - numeric_traits<T>::one();
288 if (rs > 0) {
289 // adding one artificial
290 set_type_for_logical(artificial, column_type::fixed);
291 (*this->m_A)(row, artificial) = numeric_traits<T>::one();
292 this->m_basis[row] = artificial;
293 this->m_costs[artificial] = numeric_traits<T>::zero();
294 artificial++;
295 } else {
296 // we can put a slack_var into the basis, and avoid adding an artificial variable
297 this->m_basis[row] = slack_var;
298 this->m_costs[slack_var] = numeric_traits<T>::zero();
299 }
300 slack_var++;
301 break;
302 case Less_or_equal:
303 // introduce a non-negative slack variable
304 set_type_for_logical(slack_var, column_type::lower_bound);
305 (*this->m_A)(row, slack_var) = numeric_traits<T>::one();
306 if (rs < 0) {
307 // adding one artificial
308 set_type_for_logical(artificial, column_type::fixed);
309 (*this->m_A)(row, artificial) = - numeric_traits<T>::one();
310 this->m_basis[row] = artificial;
311 this->m_costs[artificial] = numeric_traits<T>::zero();
312 artificial++;
313 } else {
314 // we can put slack_var into the basis, and avoid adding an artificial variable
315 this->m_basis[row] = slack_var;
316 this->m_costs[slack_var] = numeric_traits<T>::zero();
317 }
318 slack_var++;
319 break;
320 }
321 }
322
augment_matrix_A_and_fill_x_and_allocate_some_fields()323 template <typename T, typename X> void lp_dual_simplex<T, X>::augment_matrix_A_and_fill_x_and_allocate_some_fields() {
324 this->count_slacks_and_artificials();
325 this->m_A->add_columns_at_the_end(this->m_slacks + this->m_artificials);
326 unsigned n = this->m_A->column_count();
327 this->m_column_types_of_core_solver.resize(n);
328 m_column_types_of_logicals.resize(this->m_slacks + this->m_artificials);
329 this->m_costs.resize(n);
330 this->m_upper_bounds.resize(n);
331 this->m_lower_bounds.resize(n);
332 m_can_enter_basis.resize(n);
333 this->m_basis.resize(this->m_A->row_count());
334 }
335
336
337
copy_m_b_aside_and_set_it_to_zeros()338 template <typename T, typename X> void lp_dual_simplex<T, X>::copy_m_b_aside_and_set_it_to_zeros() {
339 for (unsigned i = 0; i < this->m_b.size(); i++) {
340 m_b_copy.push_back(this->m_b[i]);
341 this->m_b[i] = numeric_traits<T>::zero(); // preparing for the first stage
342 }
343 }
344
find_maximal_solution()345 template <typename T, typename X> void lp_dual_simplex<T, X>::find_maximal_solution(){
346 if (this->problem_is_empty()) {
347 this->m_status = lp_status::EMPTY;
348 return;
349 }
350
351 this->flip_costs(); // do it for now, todo ( remove the flipping)
352
353 this->cleanup();
354 if (this->m_status == lp_status::INFEASIBLE) {
355 return;
356 }
357 this->fill_matrix_A_and_init_right_side();
358 this->fill_m_b();
359 this->scale();
360 augment_matrix_A_and_fill_x_and_allocate_some_fields();
361 fill_first_stage_solver_fields();
362 copy_m_b_aside_and_set_it_to_zeros();
363 stage1();
364 if (this->m_status == lp_status::FEASIBLE) {
365 stage2();
366 }
367 }
368
369
get_current_cost()370 template <typename T, typename X> T lp_dual_simplex<T, X>::get_current_cost() const {
371 T ret = numeric_traits<T>::zero();
372 for (auto it : this->m_map_from_var_index_to_column_info) {
373 ret += this->get_column_cost_value(it.first, it.second);
374 }
375 return -ret; // we flip costs for now
376 }
377 }
378