25#include "absl/flags/flag.h"
26#include "absl/log/check.h"
27#include "absl/log/vlog_is_on.h"
28#include "absl/random/bit_gen_ref.h"
29#include "absl/random/random.h"
30#include "absl/random/seed_sequences.h"
31#include "absl/strings/match.h"
32#include "absl/strings/str_cat.h"
33#include "absl/strings/str_format.h"
34#include "absl/strings/string_view.h"
57ABSL_FLAG(
bool, simplex_display_numbers_as_fractions,
false,
58 "Display numbers as fractions.");
59ABSL_FLAG(
bool, simplex_stop_after_first_basis,
false,
60 "Stop after first basis has been computed.");
61ABSL_FLAG(
bool, simplex_stop_after_feasibility,
false,
62 "Stop after first phase has been completed.");
63ABSL_FLAG(
bool, simplex_display_stats,
false,
"Display algorithm statistics.");
73 explicit Cleanup(std::function<
void()> closure)
74 : closure_(
std::move(closure)) {}
75 ~Cleanup() { closure_(); }
78 std::function<void()> closure_;
82#define DCHECK_COL_BOUNDS(col) \
85 DCHECK_GT(num_cols_, col); \
89#define DCHECK_ROW_BOUNDS(row) \
92 DCHECK_GT(num_rows_, row); \
112 random_(
absl::BitGenRef(deterministic_random_)),
113 basis_factorization_(&compact_matrix_, &basis_),
114 variables_info_(compact_matrix_),
115 primal_edge_norms_(compact_matrix_, variables_info_,
116 basis_factorization_),
117 dual_edge_norms_(basis_factorization_),
118 dual_prices_(random_),
119 variable_values_(parameters_, compact_matrix_, basis_, variables_info_,
120 basis_factorization_, &dual_edge_norms_, &dual_prices_),
121 update_row_(compact_matrix_, transposed_matrix_, variables_info_, basis_,
122 basis_factorization_),
123 reduced_costs_(compact_matrix_, objective_, basis_, variables_info_,
124 basis_factorization_, random_),
125 entering_variable_(variables_info_, random_, &reduced_costs_),
126 primal_prices_(random_, variables_info_, &primal_edge_norms_,
130 function_stats_(
"SimplexFunctionStats"),
138 solution_state_.statuses.clear();
139 variable_starting_values_.clear();
147 if (state.
statuses == solution_state_.statuses)
return;
149 solution_state_ = state;
150 solution_state_has_been_set_externally_ =
true;
155 variable_starting_values_ = values;
161 const double start_time =
time_limit->GetElapsedTime();
162 default_logger_.EnableLogging(parameters_.log_search_progress());
163 default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
164 parameters_ = initial_parameters_;
165 PropagateParameters();
168 if (transpose_was_changed_) {
169 compact_matrix_.PopulateFromTranspose(transposed_matrix_);
170 num_rows_ = compact_matrix_.num_rows();
171 num_cols_ = compact_matrix_.num_cols();
175 DCHECK_EQ(num_cols_, objective.
size());
178 objective_scaling_factor_ = objective_scaling_factor;
179 objective_offset_ = objective_offset;
180 const bool objective_is_unchanged = objective_ == objective;
181 objective_ = objective;
182 InitializeObjectiveLimit();
185 variables_info_.InitializeFromMutatedState();
187 if (objective_is_unchanged && parameters_.use_dual_simplex() &&
188 !transpose_was_changed_ && !solution_state_has_been_set_externally_ &&
189 !solution_state_.IsEmpty()) {
191 primal_edge_norms_.Clear();
192 variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
194 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
195 variable_values_.RecomputeBasicVariableValues();
196 return SolveInternal(start_time,
false, objective,
time_limit);
201 return SolveInternal(start_time,
false, objective,
time_limit);
205 const double start_time =
time_limit->GetElapsedTime();
206 default_logger_.EnableLogging(parameters_.log_search_progress());
207 default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
215ABSL_MUST_USE_RESULT
Status RevisedSimplex::SolveInternal(
216 double start_time,
bool is_maximization_problem,
220 Cleanup update_deterministic_time_on_return(
228 DisplayBasicVariableStatistics();
231 dual_infeasibility_improvement_direction_.clear();
235 phase_ = Phase::FEASIBILITY;
237 num_feasibility_iterations_ = 0;
238 num_optimization_iterations_ = 0;
239 num_push_iterations_ = 0;
240 feasibility_time_ = 0.0;
241 optimization_time_ = 0.0;
249 solution_state_has_been_set_externally_ =
true;
252 ComputeNumberOfEmptyRows();
253 ComputeNumberOfEmptyColumns();
256 if (absl::GetFlag(FLAGS_simplex_stop_after_first_basis)) {
261 const bool use_dual = parameters_.use_dual_simplex();
266 primal_edge_norms_.SetPricingRule(parameters_.feasibility_rule());
268 if (parameters_.perturb_costs_in_dual_simplex()) {
269 reduced_costs_.PerturbCosts();
272 if (parameters_.use_dedicated_dual_feasibility_algorithm()) {
273 variables_info_.MakeBoxedVariableRelevant(
false);
275 DualMinimize(phase_ == Phase::FEASIBILITY,
time_limit));
288 variables_info_.MakeBoxedVariableRelevant(
true);
289 reduced_costs_.MakeReducedCostsPrecise();
292 MakeBoxedVariableDualFeasible(
293 variables_info_.GetNonBasicBoxedVariables(),
295 variable_values_.RecomputeBasicVariableValues();
301 reduced_costs_.MakeReducedCostsPrecise();
302 bool refactorize = reduced_costs_.NeedsBasisRefactorization();
306 reduced_costs_.ComputeMaximumDualInfeasibilityOnNonBoxedVariables();
307 if (initial_infeasibility <
308 reduced_costs_.GetDualFeasibilityTolerance()) {
309 SOLVER_LOG(logger_,
"Initial basis is dual feasible.");
311 MakeBoxedVariableDualFeasible(
312 variables_info_.GetNonBasicBoxedVariables(),
314 variable_values_.RecomputeBasicVariableValues();
317 variables_info_.TransformToDualPhaseIProblem(
318 reduced_costs_.GetDualFeasibilityTolerance(),
319 reduced_costs_.GetReducedCosts());
321 variable_values_.ResetAllNonBasicVariableValues(zero);
322 variable_values_.RecomputeBasicVariableValues();
330 variables_info_.EndDualPhaseI(
331 reduced_costs_.GetDualFeasibilityTolerance(),
332 reduced_costs_.GetFullReducedCosts());
333 variable_values_.ResetAllNonBasicVariableValues(
334 variable_starting_values_);
335 variable_values_.RecomputeBasicVariableValues();
345 if (reduced_costs_.ComputeMaximumDualInfeasibility() <
346 reduced_costs_.GetDualFeasibilityTolerance() + 1e-6) {
349 SOLVER_LOG(logger_,
"Infeasible after first phase.");
360 objective_ = objective_coefficients;
361 if (is_maximization_problem) {
366 objective_.
resize(num_cols_, 0.0);
367 reduced_costs_.ResetForNewObjective();
373 phase_ = Phase::OPTIMIZATION;
374 feasibility_time_ =
time_limit->GetElapsedTime() - start_time;
375 primal_edge_norms_.SetPricingRule(parameters_.optimization_rule());
376 num_feasibility_iterations_ = num_iterations_;
391 for (
int num_optims = 0;
395 num_optims <= parameters_.max_number_of_reoptimizations() &&
396 !objective_limit_reached_ &&
397 (num_iterations_ == 0 ||
398 num_iterations_ < parameters_.max_number_of_iterations()) &&
400 !absl::GetFlag(FLAGS_simplex_stop_after_feasibility) &&
410 DualMinimize(phase_ == Phase::FEASIBILITY,
time_limit));
418 basis_factorization_.IsRefactorized());
421 if (!integrality_scale_.empty() &&
430 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
433 variable_values_.RecomputeBasicVariableValues();
434 reduced_costs_.ClearAndRemoveCostShifts();
448 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
449 if (reduced_costs_.ComputeMaximumDualResidual() > tolerance ||
450 variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
451 variable_values_.ComputeMaximumPrimalInfeasibility() > tolerance) {
453 "PRIMAL_UNBOUNDED was reported, but the residual and/or "
454 "dual infeasibility is above the tolerance");
455 if (parameters_.change_status_to_imprecise()) {
478 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
479 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
481 for (ColIndex col(0); col < num_cols_; ++col) {
482 cost_delta += solution_primal_ray_[col] * objective_[col];
483 if (solution_primal_ray_[col] > 0 && upper_bounds[col] !=
kInfinity) {
484 const Fractional value = variable_values_.Get(col);
485 const Fractional distance = (upper_bounds[col] - value + tolerance) /
486 solution_primal_ray_[col];
487 min_distance = std::min(distance, min_distance);
488 max_magnitude = std::max(solution_primal_ray_[col], max_magnitude);
490 if (solution_primal_ray_[col] < 0 && lower_bounds[col] != -
kInfinity) {
491 const Fractional value = variable_values_.Get(col);
492 const Fractional distance = (value - lower_bounds[col] + tolerance) /
493 -solution_primal_ray_[col];
494 min_distance = std::min(distance, min_distance);
495 max_magnitude = std::max(-solution_primal_ray_[col], max_magnitude);
498 SOLVER_LOG(logger_,
"Primal unbounded ray: max blocking magnitude = ",
499 max_magnitude,
", min distance to bound + ", tolerance,
" = ",
500 min_distance,
", ray cost delta = ", cost_delta);
501 if (min_distance * std::abs(cost_delta) < 1 &&
502 reduced_costs_.ComputeMaximumDualInfeasibility() <= tolerance) {
504 "PRIMAL_UNBOUNDED was reported, but the tolerance are good "
505 "and the unbounded ray is not great.");
507 "The difference between unbounded and optimal can depend "
508 "on a slight change of tolerance, trying to see if we are "
509 "at OPTIMAL after postsolve.");
515 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
516 if (reduced_costs_.ComputeMaximumDualResidual() > tolerance ||
517 variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
518 reduced_costs_.ComputeMaximumDualInfeasibility() > tolerance) {
520 "DUAL_UNBOUNDED was reported, but the residual and/or "
521 "dual infeasibility is above the tolerance");
522 if (parameters_.change_status_to_imprecise()) {
531 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
532 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
535 for (ColIndex col(0); col < num_cols_; ++col) {
536 const Fractional coeff = solution_dual_ray_row_combination_[col];
539 error = std::max(error, coeff);
541 implied_lb += coeff * lower_bounds[col];
543 }
else if (coeff < 0) {
545 error = std::max(error, -coeff);
547 implied_lb += coeff * upper_bounds[col];
552 " infeasibility=", implied_lb);
553 if (implied_lb < tolerance || error > tolerance) {
555 "DUAL_UNBOUNDED was reported, but the dual ray is not "
556 "proving infeasibility with high enough tolerance");
557 if (parameters_.change_status_to_imprecise()) {
568 parameters_.solution_feasibility_tolerance();
570 variable_values_.ComputeMaximumPrimalResidual();
572 reduced_costs_.ComputeMaximumDualResidual();
573 if (primal_residual > solution_tolerance ||
574 dual_residual > solution_tolerance) {
576 "OPTIMAL was reported, yet one of the residuals is "
577 "above the solution feasibility tolerance after the "
578 "shift/perturbation are removed.");
579 if (parameters_.change_status_to_imprecise()) {
589 std::max(primal_residual,
590 Fractional(parameters_.primal_feasibility_tolerance()));
592 std::max(dual_residual,
593 Fractional(parameters_.dual_feasibility_tolerance()));
595 variable_values_.ComputeMaximumPrimalInfeasibility();
597 reduced_costs_.ComputeMaximumDualInfeasibility();
598 if (primal_infeasibility > primal_tolerance &&
599 dual_infeasibility > dual_tolerance) {
601 "OPTIMAL was reported, yet both of the infeasibility "
602 "are above the tolerance after the "
603 "shift/perturbation are removed.");
604 if (parameters_.change_status_to_imprecise()) {
607 }
else if (primal_infeasibility > primal_tolerance) {
608 if (num_optims == parameters_.max_number_of_reoptimizations()) {
610 "The primal infeasibility is still higher than the "
611 "requested internal tolerance, but the maximum "
612 "number of optimization is reached.");
616 SOLVER_LOG(logger_,
"Re-optimizing with dual simplex ... ");
618 }
else if (dual_infeasibility > dual_tolerance) {
619 if (num_optims == parameters_.max_number_of_reoptimizations()) {
621 "The dual infeasibility is still higher than the "
622 "requested internal tolerance, but the maximum "
623 "number of optimization is reached.");
627 SOLVER_LOG(logger_,
"Re-optimizing with primal simplex ... ");
638 if (parameters_.change_status_to_imprecise() &&
640 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
641 if (variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
642 reduced_costs_.ComputeMaximumDualResidual() > tolerance) {
647 if (reduced_costs_.ComputeMaximumDualInfeasibility() > tolerance) {
653 if (variable_values_.ComputeMaximumPrimalInfeasibility() > tolerance) {
659 total_time_ =
time_limit->GetElapsedTime() - start_time;
660 optimization_time_ = total_time_ - feasibility_time_;
661 num_optimization_iterations_ = num_iterations_ - num_feasibility_iterations_;
665 if (!variable_starting_values_.empty()) {
666 const int num_super_basic = ComputeNumberOfSuperBasicVariables();
667 if (num_super_basic > 0) {
669 "Num super-basic variables left after optimize phase: ",
671 if (parameters_.push_to_vertex()) {
674 phase_ = Phase::PUSH;
680 "Skipping push phase because optimize didn't succeed.");
686 total_time_ =
time_limit->GetElapsedTime() - start_time;
687 push_time_ = total_time_ - feasibility_time_ - optimization_time_;
688 num_push_iterations_ = num_iterations_ - num_feasibility_iterations_ -
689 num_optimization_iterations_;
692 solution_objective_value_ = ComputeInitialProblemObjectiveValue();
693 solution_dual_values_ = reduced_costs_.GetDualValues();
694 solution_reduced_costs_ = reduced_costs_.GetReducedCosts();
697 if (is_maximization_problem) {
705 solution_objective_value_ =
708 if (is_maximization_problem) {
709 solution_objective_value_ = -solution_objective_value_;
713 variable_starting_values_.clear();
719 return problem_status_;
723 return solution_objective_value_;
727 return num_iterations_;
735 return variable_values_.Get(col);
739 return solution_reduced_costs_[col];
743 return solution_reduced_costs_;
747 return solution_dual_values_[row];
751 return variables_info_.GetStatusRow()[col];
759 return -variable_values_.Get(SlackColIndex(row));
765 const VariableStatus s = variables_info_.GetStatusRow()[SlackColIndex(row)];
777 return solution_primal_ray_;
781 return solution_dual_ray_;
786 return solution_dual_ray_row_combination_;
792 DCHECK(basis_factorization_.GetColumnPermutation().empty());
793 return basis_factorization_;
796std::string RevisedSimplex::GetPrettySolverStats()
const {
797 return absl::StrFormat(
798 "Problem status : %s\n"
799 "Solving time : %-6.4g\n"
800 "Number of iterations : %u\n"
801 "Time for solvability (first phase) : %-6.4g\n"
802 "Number of iterations for solvability : %u\n"
803 "Time for optimization : %-6.4g\n"
804 "Number of iterations for optimization : %u\n"
805 "Stop after first basis : %d\n",
807 feasibility_time_, num_feasibility_iterations_, optimization_time_,
808 num_optimization_iterations_,
809 absl::GetFlag(FLAGS_simplex_stop_after_first_basis));
815 basis_factorization_.DeterministicTime() +
816 update_row_.DeterministicTime() +
817 entering_variable_.DeterministicTime() +
818 reduced_costs_.DeterministicTime() +
819 primal_edge_norms_.DeterministicTime();
822void RevisedSimplex::SetVariableNames() {
823 variable_name_.
resize(num_cols_,
"");
824 for (ColIndex col(0); col < first_slack_col_; ++col) {
825 const ColIndex var_index = col + 1;
826 variable_name_[col] = absl::StrFormat(
"x%d",
ColToIntIndex(var_index));
828 for (ColIndex col(first_slack_col_); col < num_cols_; ++col) {
829 const ColIndex var_index = col - first_slack_col_ + 1;
830 variable_name_[col] = absl::StrFormat(
"s%d",
ColToIntIndex(var_index));
834void RevisedSimplex::SetNonBasicVariableStatusAndDeriveValue(
836 variables_info_.UpdateToNonBasicStatus(col, status);
837 variable_values_.SetNonBasicVariableValueFromStatus(col);
840bool RevisedSimplex::BasisIsConsistent()
const {
841 const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow();
843 for (RowIndex row(0); row < num_rows_; ++row) {
844 const ColIndex col = basis_[row];
845 if (!is_basic.IsSet(col))
return false;
848 ColIndex cols_in_basis(0);
849 ColIndex cols_not_in_basis(0);
850 for (ColIndex col(0); col < num_cols_; ++col) {
851 cols_in_basis += is_basic.IsSet(col);
852 cols_not_in_basis += !is_basic.IsSet(col);
853 if (is_basic.IsSet(col) !=
859 if (cols_not_in_basis != num_cols_ -
RowToColIndex(num_rows_))
return false;
865void RevisedSimplex::UpdateBasis(ColIndex entering_col, RowIndex basis_row,
873 DCHECK(!variables_info_.GetIsBasicBitRow().IsSet(entering_col));
874 DCHECK_NE(basis_[basis_row], entering_col);
877 const ColIndex leaving_col = basis_[basis_row];
878 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(leaving_col));
883 variables_info_.UpdateToNonBasicStatus(leaving_col, leaving_variable_status);
888 basis_[basis_row] = entering_col;
889 variables_info_.UpdateToBasicStatus(entering_col);
890 update_row_.Invalidate();
896class ColumnComparator {
898 explicit ColumnComparator(
const DenseRow& value) : value_(value) {}
899 bool operator()(ColIndex col_a, ColIndex col_b)
const {
900 return value_[col_a] < value_[col_b];
917void RevisedSimplex::UseSingletonColumnInInitialBasis(
RowToColMapping* basis) {
924 std::vector<ColIndex> singleton_column;
925 DenseRow cost_variation(num_cols_, 0.0);
926 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
927 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
928 for (ColIndex col(0); col < num_cols_; ++col) {
929 if (compact_matrix_.column(col).num_entries() != 1)
continue;
930 if (lower_bounds[col] == upper_bounds[col])
continue;
931 const Fractional slope = compact_matrix_.column(col).GetFirstCoefficient();
932 if (variable_values_.Get(col) == lower_bounds[col]) {
933 cost_variation[col] = objective_[col] / std::abs(slope);
935 cost_variation[col] = -objective_[col] / std::abs(slope);
937 singleton_column.push_back(col);
939 if (singleton_column.empty())
return;
946 ColumnComparator comparator(cost_variation);
947 std::sort(singleton_column.begin(), singleton_column.end(), comparator);
948 DCHECK_LE(cost_variation[singleton_column.front()],
949 cost_variation[singleton_column.back()]);
956 const DenseRow& variable_values = variable_values_.GetDenseRow();
957 for (
const ColIndex col : singleton_column) {
958 const RowIndex row = compact_matrix_.column(col).EntryRow(EntryIndex(0));
969 if (error_[row] == 0.0)
continue;
975 compact_matrix_.column(col).EntryCoefficient(EntryIndex(0));
976 const Fractional new_value = variable_values[col] + error_[row] / coeff;
977 if (new_value >= lower_bounds[col] && new_value <= upper_bounds[col]) {
988 const Fractional box_width = variables_info_.GetBoundDifference(col);
989 DCHECK_NE(box_width, 0.0);
990 DCHECK_NE(error_[row], 0.0);
991 const Fractional error_sign = error_[row] / coeff;
992 if (variable_values[col] == lower_bounds[col] && error_sign > 0.0) {
994 error_[row] -= coeff * box_width;
995 SetNonBasicVariableStatusAndDeriveValue(col,
999 if (variable_values[col] == upper_bounds[col] && error_sign < 0.0) {
1001 error_[row] += coeff * box_width;
1002 SetNonBasicVariableStatusAndDeriveValue(col,
1009bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged(
1011 bool* only_change_is_new_rows,
bool* only_change_is_new_cols,
1012 ColIndex* num_new_cols) {
1014 DCHECK(only_change_is_new_rows !=
nullptr);
1015 DCHECK(only_change_is_new_cols !=
nullptr);
1016 DCHECK(num_new_cols !=
nullptr);
1017 DCHECK_EQ(num_cols_, compact_matrix_.num_cols());
1018 DCHECK_EQ(num_rows_, compact_matrix_.num_rows());
1021 const bool old_part_of_matrix_is_unchanged =
1023 num_rows_, first_slack_col_, lp.GetSparseMatrix(), compact_matrix_);
1026 const ColIndex lp_first_slack =
1027 lp_is_in_equation_form ? lp.GetFirstSlackVariable() : lp.num_variables();
1032 if (old_part_of_matrix_is_unchanged && lp.num_constraints() == num_rows_ &&
1033 lp_first_slack == first_slack_col_) {
1038 if (parameters_.use_transposed_matrix()) {
1039 if (transposed_matrix_.IsEmpty()) {
1040 transposed_matrix_.PopulateFromTranspose(compact_matrix_);
1043 transposed_matrix_.Reset(RowIndex(0));
1050 *only_change_is_new_rows = old_part_of_matrix_is_unchanged &&
1051 lp.num_constraints() > num_rows_ &&
1052 lp_first_slack == first_slack_col_;
1056 *only_change_is_new_cols = old_part_of_matrix_is_unchanged &&
1057 lp.num_constraints() == num_rows_ &&
1058 lp_first_slack > first_slack_col_;
1059 *num_new_cols = *only_change_is_new_cols ? lp_first_slack - first_slack_col_
1063 first_slack_col_ = lp_first_slack;
1066 num_rows_ = lp.num_constraints();
1067 num_cols_ = lp_first_slack +
RowToColIndex(lp.num_constraints());
1070 if (lp_is_in_equation_form) {
1073 compact_matrix_.PopulateFromMatrixView(MatrixView(lp.GetSparseMatrix()));
1075 compact_matrix_.PopulateFromSparseMatrixAndAddSlacks(lp.GetSparseMatrix());
1077 if (parameters_.use_transposed_matrix()) {
1078 transposed_matrix_.PopulateFromTranspose(compact_matrix_);
1080 transposed_matrix_.Reset(RowIndex(0));
1087bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
1089 ColIndex num_new_cols) {
1091 DCHECK_LE(num_new_cols, first_slack_col_);
1092 const ColIndex first_new_col(first_slack_col_ - num_new_cols);
1095 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1096 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1097 for (ColIndex col(0); col < first_new_col; ++col) {
1098 if (lower_bounds[col] != lp.variable_lower_bounds()[col] ||
1099 upper_bounds[col] != lp.variable_upper_bounds()[col]) {
1105 for (ColIndex col(first_new_col); col < first_slack_col_; ++col) {
1106 if (lp.variable_lower_bounds()[col] != 0.0 &&
1107 lp.variable_upper_bounds()[col] != 0.0) {
1113 if (lp_is_in_equation_form) {
1114 for (ColIndex col(first_slack_col_); col < num_cols_; ++col) {
1115 if (lower_bounds[col - num_new_cols] != lp.variable_lower_bounds()[col] ||
1116 upper_bounds[col - num_new_cols] != lp.variable_upper_bounds()[col]) {
1121 DCHECK_EQ(num_rows_, lp.num_constraints());
1122 for (RowIndex row(0); row < num_rows_; ++row) {
1124 if (lower_bounds[col - num_new_cols] !=
1125 -lp.constraint_upper_bounds()[row] ||
1126 upper_bounds[col - num_new_cols] !=
1127 -lp.constraint_lower_bounds()[row]) {
1135bool RevisedSimplex::InitializeObjectiveAndTestIfUnchanged(
1139 bool objective_is_unchanged =
true;
1140 objective_.resize(num_cols_, 0.0);
1144 DCHECK_GE(num_cols_, lp.num_variables());
1146 const auto obj_coeffs = lp.objective_coefficients().const_view();
1147 for (ColIndex col(obj_coeffs.size()); col < num_cols_; ++col) {
1148 if (objective_[col] != 0.0) {
1149 objective_is_unchanged =
false;
1150 objective_[col] = 0.0;
1154 if (lp.IsMaximizationProblem()) {
1156 for (ColIndex col(0); col < obj_coeffs.size(); ++col) {
1158 if (objective_[col] != coeff) {
1159 objective_is_unchanged =
false;
1160 objective_[col] = coeff;
1163 objective_offset_ = -lp.objective_offset();
1164 objective_scaling_factor_ = -lp.objective_scaling_factor();
1166 for (ColIndex col(0); col < obj_coeffs.size(); ++col) {
1168 if (objective_[col] != coeff) {
1169 objective_is_unchanged =
false;
1170 objective_[col] = coeff;
1173 objective_offset_ = lp.objective_offset();
1174 objective_scaling_factor_ = lp.objective_scaling_factor();
1177 return objective_is_unchanged;
1180void RevisedSimplex::InitializeObjectiveLimit() {
1181 objective_limit_reached_ =
false;
1182 DCHECK(std::isfinite(objective_offset_));
1183 DCHECK(std::isfinite(objective_scaling_factor_));
1184 DCHECK_NE(0.0, objective_scaling_factor_);
1187 for (
const bool set_dual : {
true,
false}) {
1199 const Fractional limit = (objective_scaling_factor_ >= 0.0) != set_dual
1200 ? parameters_.objective_lower_limit()
1201 : parameters_.objective_upper_limit();
1203 limit / objective_scaling_factor_ - objective_offset_;
1205 dual_objective_limit_ = shifted_limit;
1207 primal_objective_limit_ = shifted_limit;
1217Status RevisedSimplex::CreateInitialBasis() {
1224 variables_info_.InitializeToDefaultStatus();
1225 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
1229 for (RowIndex row(0); row < num_rows_; ++row) {
1230 basis[row] = SlackColIndex(row);
1235 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1236 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1237 if (!parameters_.use_dual_simplex() &&
1239 parameters_.exploit_singleton_column_in_initial_basis()) {
1243 for (ColIndex col(0); col < num_cols_; ++col) {
1244 if (compact_matrix_.column(col).num_entries() != 1)
continue;
1245 const VariableStatus status = variables_info_.GetStatusRow()[col];
1246 const Fractional objective = objective_[col];
1247 if (objective > 0 &&
IsFinite(lower_bounds[col]) &&
1249 SetNonBasicVariableStatusAndDeriveValue(col,
1251 }
else if (objective < 0 &&
IsFinite(upper_bounds[col]) &&
1253 SetNonBasicVariableStatusAndDeriveValue(col,
1260 ComputeVariableValuesError();
1269 UseSingletonColumnInInitialBasis(&basis);
1272 for (RowIndex row(0); row < num_rows_; ++row) {
1274 basis[row] = SlackColIndex(row);
1281 return InitializeFirstBasis(basis);
1284 InitialBasis initial_basis(compact_matrix_, objective_, lower_bounds,
1285 upper_bounds, variables_info_.GetTypeRow());
1286 if (parameters_.use_dual_simplex()) {
1289 initial_basis.GetDualMarosBasis(num_cols_, &basis);
1291 initial_basis.GetPrimalMarosBasis(num_cols_, &basis);
1293 int number_changed = 0;
1294 for (RowIndex row(0); row < num_rows_; ++row) {
1295 if (basis[row] != SlackColIndex(row)) {
1299 VLOG(1) <<
"Number of Maros basis changes: " << number_changed;
1303 int num_fixed_variables = 0;
1304 for (RowIndex row(0); row < basis.size(); ++row) {
1305 const ColIndex col = basis[row];
1306 if (lower_bounds[col] == upper_bounds[col]) {
1308 ++num_fixed_variables;
1312 if (num_fixed_variables == 0) {
1313 SOLVER_LOG(logger_,
"Crash is set to ", parameters_.initial_basis(),
1314 " but there is no equality rows to remove from initial all "
1315 "slack basis. Starting from there.");
1318 SOLVER_LOG(logger_,
"Trying to remove ", num_fixed_variables,
1319 " fixed variables from the initial basis.");
1320 InitialBasis initial_basis(compact_matrix_, objective_, lower_bounds,
1321 upper_bounds, variables_info_.GetTypeRow());
1324 if (parameters_.use_scaling()) {
1325 initial_basis.CompleteBixbyBasis(first_slack_col_, &basis);
1327 VLOG(1) <<
"Bixby initial basis algorithm requires the problem "
1328 <<
"to be scaled. Skipping Bixby's algorithm.";
1333 if (parameters_.use_dual_simplex()) {
1336 initial_basis.CompleteTriangularDualBasis(num_cols_, &basis);
1338 initial_basis.CompleteTriangularPrimalBasis(num_cols_, &basis);
1341 const Status status = InitializeFirstBasis(basis);
1347 "Advanced basis algo failed, Reverting to all slack basis.");
1349 for (RowIndex row(0); row < num_rows_; ++row) {
1350 basis[row] = SlackColIndex(row);
1356 LOG(WARNING) <<
"Unsupported initial_basis parameters: "
1357 << parameters_.initial_basis();
1360 return InitializeFirstBasis(basis);
1369 for (RowIndex row(0); row < num_rows_; ++row) {
1371 basis_[row] = SlackColIndex(row);
1384 basis_factorization_.ComputeInfinityNormConditionNumberUpperBound();
1385 if (condition_number_ub > parameters_.initial_condition_number_threshold()) {
1386 const std::string error_message =
1387 absl::StrCat(
"The matrix condition number upper bound is too high: ",
1388 condition_number_ub);
1394 for (RowIndex row(0); row < num_rows_; ++row) {
1395 variables_info_.UpdateToBasicStatus(basis_[row]);
1397 DCHECK(BasisIsConsistent());
1399 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
1400 variable_values_.RecomputeBasicVariableValues();
1402 if (logger_->LoggingIsEnabled()) {
1406 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
1407 if (variable_values_.ComputeMaximumPrimalResidual() > tolerance) {
1410 "The primal residual of the initial basis is above the tolerance, ",
1411 variable_values_.ComputeMaximumPrimalResidual(),
" vs. ", tolerance);
1418 parameters_ = initial_parameters_;
1419 PropagateParameters();
1427 const bool lp_is_in_equation_form = lp.IsInEquationForm();
1434 ColIndex num_new_cols(0);
1435 bool only_change_is_new_rows =
false;
1436 bool only_change_is_new_cols =
false;
1437 const bool matrix_is_unchanged = InitializeMatrixAndTestIfUnchanged(
1438 lp, lp_is_in_equation_form, &only_change_is_new_rows,
1439 &only_change_is_new_cols, &num_new_cols);
1440 const bool only_new_bounds =
1441 only_change_is_new_cols && num_new_cols > 0 &&
1442 OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
1443 lp, lp_is_in_equation_form, num_new_cols);
1446 const bool objective_is_unchanged = InitializeObjectiveAndTestIfUnchanged(lp);
1448 const bool bounds_are_unchanged =
1449 lp_is_in_equation_form
1450 ? variables_info_.LoadBoundsAndReturnTrueIfUnchanged(
1451 lp.variable_lower_bounds(), lp.variable_upper_bounds())
1452 : variables_info_.LoadBoundsAndReturnTrueIfUnchanged(
1453 lp.variable_lower_bounds(), lp.variable_upper_bounds(),
1454 lp.constraint_lower_bounds(), lp.constraint_upper_bounds());
1459 if (matrix_is_unchanged && parameters_.allow_simplex_algorithm_change()) {
1460 if (objective_is_unchanged && !bounds_are_unchanged) {
1461 parameters_.set_use_dual_simplex(
true);
1462 PropagateParameters();
1464 if (bounds_are_unchanged && !objective_is_unchanged) {
1465 parameters_.set_use_dual_simplex(
false);
1466 PropagateParameters();
1470 InitializeObjectiveLimit();
1475 if (VLOG_IS_ON(2)) {
1486 bool solve_from_scratch =
true;
1489 if (!solution_state_.IsEmpty() && !solution_state_has_been_set_externally_) {
1490 if (!parameters_.use_dual_simplex()) {
1495 dual_edge_norms_.Clear();
1496 dual_pricing_vector_.clear();
1497 if (matrix_is_unchanged && bounds_are_unchanged) {
1500 reduced_costs_.ClearAndRemoveCostShifts();
1501 solve_from_scratch =
false;
1502 }
else if (only_change_is_new_cols && only_new_bounds) {
1503 variables_info_.InitializeFromBasisState(first_slack_col_, num_new_cols,
1505 variable_values_.ResetAllNonBasicVariableValues(
1506 variable_starting_values_);
1508 const ColIndex first_new_col(first_slack_col_ - num_new_cols);
1509 for (ColIndex& col_ref : basis_) {
1510 if (col_ref >= first_new_col) {
1511 col_ref += num_new_cols;
1518 primal_edge_norms_.Clear();
1519 reduced_costs_.ClearAndRemoveCostShifts();
1520 solve_from_scratch =
false;
1526 primal_edge_norms_.Clear();
1527 if (objective_is_unchanged) {
1528 if (matrix_is_unchanged) {
1529 if (!bounds_are_unchanged) {
1530 variables_info_.InitializeFromBasisState(
1531 first_slack_col_, ColIndex(0), solution_state_);
1532 variable_values_.ResetAllNonBasicVariableValues(
1533 variable_starting_values_);
1534 variable_values_.RecomputeBasicVariableValues();
1536 solve_from_scratch =
false;
1537 }
else if (only_change_is_new_rows) {
1540 variables_info_.InitializeFromBasisState(
1541 first_slack_col_, ColIndex(0), solution_state_);
1542 dual_edge_norms_.ResizeOnNewRows(num_rows_);
1547 reduced_costs_.ClearAndRemoveCostShifts();
1548 dual_pricing_vector_.clear();
1551 if (InitializeFirstBasis(basis_).ok()) {
1552 solve_from_scratch =
false;
1559 return FinishInitialization(solve_from_scratch);
1562Status RevisedSimplex::FinishInitialization(
bool solve_from_scratch) {
1565 if (solve_from_scratch && !solution_state_.IsEmpty()) {
1566 basis_factorization_.Clear();
1567 reduced_costs_.ClearAndRemoveCostShifts();
1568 primal_edge_norms_.Clear();
1569 dual_edge_norms_.Clear();
1570 dual_pricing_vector_.clear();
1574 variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
1578 std::vector<ColIndex> candidates;
1579 for (
const ColIndex col : variables_info_.GetIsBasicBitRow()) {
1580 candidates.push_back(col);
1582 SOLVER_LOG(logger_,
"The warm-start state contains ", candidates.size(),
1583 " candidates for the basis (num_rows = ", num_rows_.value(),
1589 if (candidates.size() == num_rows_) {
1591 for (
const ColIndex col : candidates) {
1592 basis_.push_back(col);
1598 if (InitializeFirstBasis(basis_).ok()) {
1599 solve_from_scratch =
false;
1603 if (solve_from_scratch) {
1604 basis_ = basis_factorization_.ComputeInitialBasis(candidates);
1605 const int num_super_basic =
1606 variables_info_.ChangeUnusedBasicVariablesToFree(basis_);
1607 const int num_snapped = variables_info_.SnapFreeVariablesToBound(
1608 parameters_.crossover_bound_snapping_distance(),
1609 variable_starting_values_);
1610 if (logger_->LoggingIsEnabled()) {
1611 SOLVER_LOG(logger_,
"The initial basis did not use ",
1612 " BASIC columns from the initial state and used ",
1613 (num_rows_ - (candidates.size() - num_super_basic)).value(),
1614 " slack variables that were not marked BASIC.");
1615 if (num_snapped > 0) {
1617 " of the FREE variables where moved to their bound.");
1621 if (InitializeFirstBasis(basis_).ok()) {
1622 solve_from_scratch =
false;
1625 "RevisedSimplex is not using the warm start "
1626 "basis because it is not factorizable.");
1631 if (solve_from_scratch) {
1632 SOLVER_LOG(logger_,
"Starting basis: create from scratch.");
1633 basis_factorization_.Clear();
1634 reduced_costs_.ClearAndRemoveCostShifts();
1635 primal_edge_norms_.Clear();
1636 dual_edge_norms_.Clear();
1637 dual_pricing_vector_.clear();
1640 SOLVER_LOG(logger_,
"Starting basis: incremental solve.");
1642 DCHECK(BasisIsConsistent());
1644 transpose_was_changed_ =
false;
1648void RevisedSimplex::DisplayBasicVariableStatistics() {
1651 int num_fixed_variables = 0;
1652 int num_free_variables = 0;
1653 int num_variables_at_bound = 0;
1654 int num_slack_variables = 0;
1655 int num_infeasible_variables = 0;
1657 const DenseRow& variable_values = variable_values_.GetDenseRow();
1659 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1660 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1661 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
1662 for (RowIndex row(0); row < num_rows_; ++row) {
1663 const ColIndex col = basis_[row];
1664 const Fractional value = variable_values[col];
1666 ++num_free_variables;
1668 if (value > upper_bounds[col] + tolerance ||
1669 value < lower_bounds[col] - tolerance) {
1670 ++num_infeasible_variables;
1672 if (col >= first_slack_col_) {
1673 ++num_slack_variables;
1675 if (lower_bounds[col] == upper_bounds[col]) {
1676 ++num_fixed_variables;
1677 }
else if (variable_values[col] == lower_bounds[col] ||
1678 variable_values[col] == upper_bounds[col]) {
1679 ++num_variables_at_bound;
1683 SOLVER_LOG(logger_,
"The matrix with slacks has ",
1684 compact_matrix_.num_rows().value(),
" rows, ",
1685 compact_matrix_.num_cols().value(),
" columns, ",
1686 compact_matrix_.num_entries().value(),
" entries.");
1687 SOLVER_LOG(logger_,
"Number of basic infeasible variables: ",
1688 num_infeasible_variables);
1689 SOLVER_LOG(logger_,
"Number of basic slack variables: ", num_slack_variables);
1691 "Number of basic variables at bound: ", num_variables_at_bound);
1692 SOLVER_LOG(logger_,
"Number of basic fixed variables: ", num_fixed_variables);
1693 SOLVER_LOG(logger_,
"Number of basic free variables: ", num_free_variables);
1694 SOLVER_LOG(logger_,
"Number of super-basic variables: ",
1695 ComputeNumberOfSuperBasicVariables());
1698void RevisedSimplex::SaveState() {
1699 DCHECK_EQ(num_cols_, variables_info_.GetStatusRow().size());
1700 solution_state_.statuses = variables_info_.GetStatusRow();
1701 solution_state_has_been_set_externally_ =
false;
1704RowIndex RevisedSimplex::ComputeNumberOfEmptyRows() {
1706 for (ColIndex col(0); col < num_cols_; ++col) {
1707 for (
const SparseColumn::Entry e : compact_matrix_.column(col)) {
1708 contains_data[e.row()] =
true;
1711 RowIndex num_empty_rows(0);
1712 for (RowIndex row(0); row < num_rows_; ++row) {
1713 if (!contains_data[row]) {
1715 VLOG(2) <<
"Row " << row <<
" is empty.";
1718 return num_empty_rows;
1721ColIndex RevisedSimplex::ComputeNumberOfEmptyColumns() {
1722 ColIndex num_empty_cols(0);
1723 for (ColIndex col(0); col < num_cols_; ++col) {
1724 if (compact_matrix_.column(col).IsEmpty()) {
1726 VLOG(2) <<
"Column " << col <<
" is empty.";
1729 return num_empty_cols;
1732int RevisedSimplex::ComputeNumberOfSuperBasicVariables()
const {
1734 int num_super_basic = 0;
1735 for (ColIndex col(0); col < num_cols_; ++col) {
1737 variable_values_.Get(col) != 0.0) {
1741 return num_super_basic;
1744void RevisedSimplex::CorrectErrorsOnVariableValues() {
1746 DCHECK(basis_factorization_.IsRefactorized());
1752 variable_values_.ComputeMaximumPrimalResidual();
1756 if (primal_residual >= parameters_.harris_tolerance_ratio() *
1757 parameters_.primal_feasibility_tolerance()) {
1758 variable_values_.RecomputeBasicVariableValues();
1759 VLOG(1) <<
"Primal infeasibility (bounds error) = "
1760 << variable_values_.ComputeMaximumPrimalInfeasibility()
1761 <<
", Primal residual |A.x - b| = "
1762 << variable_values_.ComputeMaximumPrimalResidual();
1766void RevisedSimplex::ComputeVariableValuesError() {
1768 error_.AssignToZero(num_rows_);
1769 const DenseRow& variable_values = variable_values_.GetDenseRow();
1770 for (ColIndex col(0); col < num_cols_; ++col) {
1771 const Fractional value = variable_values[col];
1772 compact_matrix_.ColumnAddMultipleToDenseColumn(col, -value, &error_);
1776void RevisedSimplex::ComputeDirection(ColIndex col) {
1779 basis_factorization_.RightSolveForProblemColumn(col, &direction_);
1781 if (direction_.non_zeros.empty()) {
1783 const RowIndex num_rows = num_rows_;
1785 for (RowIndex row(0); row < num_rows; ++row) {
1788 direction_.non_zeros.push_back(row);
1789 norm = std::max(norm, std::abs(value));
1793 for (
const auto e : direction_) {
1794 norm = std::max(norm, std::abs(e.coefficient()));
1797 direction_infinity_norm_ = norm;
1799 num_rows_ == 0 ? 0.0
1800 :
static_cast<double>(direction_.non_zeros.size()) /
1801 static_cast<double>(num_rows_.value())));
1804Fractional RevisedSimplex::ComputeDirectionError(ColIndex col) {
1806 compact_matrix_.ColumnCopyToDenseColumn(col, &error_);
1807 for (
const auto e : direction_) {
1808 compact_matrix_.ColumnAddMultipleToDenseColumn(col, -e.coefficient(),
1814template <
bool is_entering_reduced_cost_positive>
1817 RowIndex row)
const {
1818 const ColIndex col = basis_[row];
1819 const Fractional direction = direction_[row];
1820 const Fractional value = variable_values_.Get(col);
1821 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(col));
1822 DCHECK_NE(direction, 0.0);
1823 if (is_entering_reduced_cost_positive) {
1824 if (direction > 0.0) {
1825 return (upper_bounds[col] - value) / direction;
1827 return (lower_bounds[col] - value) / direction;
1830 if (direction > 0.0) {
1831 return (value - lower_bounds[col]) / direction;
1833 return (value - upper_bounds[col]) / direction;
1838template <
bool is_entering_reduced_cost_positive>
1839Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates(
1843 parameters_.harris_tolerance_ratio() *
1844 parameters_.primal_feasibility_tolerance();
1845 const Fractional minimum_delta = parameters_.degenerate_ministep_factor() *
1846 parameters_.primal_feasibility_tolerance();
1852 leaving_candidates->Clear();
1858 const Fractional threshold = basis_factorization_.IsRefactorized()
1859 ? parameters_.minimum_acceptable_pivot()
1860 : parameters_.ratio_test_zero_threshold();
1862 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1863 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1864 for (
const auto e : direction_) {
1865 const Fractional magnitude = std::abs(e.coefficient());
1866 if (magnitude <= threshold)
continue;
1867 const Fractional ratio = GetRatio<is_entering_reduced_cost_positive>(
1868 lower_bounds, upper_bounds, e.row());
1869 if (ratio <= harris_ratio) {
1870 leaving_candidates->SetCoefficient(e.row(), ratio);
1882 harris_ratio = std::min(harris_ratio,
1883 std::max(minimum_delta / magnitude,
1884 ratio + harris_tolerance / magnitude));
1887 return harris_ratio;
1900 if (current >= 0.0) {
1901 return candidate >= 0.0 && candidate <= current;
1903 return candidate >= current;
1911Status RevisedSimplex::ChooseLeavingVariableRow(
1912 ColIndex entering_col,
Fractional reduced_cost,
bool* refactorize,
1919 DCHECK_NE(0.0, reduced_cost);
1922 int stats_num_leaving_choices = 0;
1923 equivalent_leaving_choices_.clear();
1924 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1925 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1927 stats_num_leaving_choices = 0;
1931 const Fractional entering_value = variable_values_.Get(entering_col);
1933 (reduced_cost > 0.0) ? entering_value - lower_bounds[entering_col]
1934 : upper_bounds[entering_col] - entering_value;
1935 DCHECK_GT(current_ratio, 0.0);
1941 (reduced_cost > 0.0) ? ComputeHarrisRatioAndLeavingCandidates<true>(
1942 current_ratio, &leaving_candidates_)
1943 : ComputeHarrisRatioAndLeavingCandidates<false>(
1944 current_ratio, &leaving_candidates_);
1949 if (current_ratio <= harris_ratio) {
1951 *step_length = current_ratio;
1961 stats_num_leaving_choices = 0;
1963 equivalent_leaving_choices_.clear();
1964 for (
const SparseColumn::Entry e : leaving_candidates_) {
1966 if (ratio > harris_ratio)
continue;
1967 ++stats_num_leaving_choices;
1968 const RowIndex row = e.row();
1973 const Fractional candidate_magnitude = std::abs(direction_[row]);
1974 if (candidate_magnitude < pivot_magnitude)
continue;
1975 if (candidate_magnitude == pivot_magnitude) {
1976 if (!IsRatioMoreOrEquallyStable(ratio, current_ratio))
continue;
1977 if (ratio == current_ratio) {
1979 equivalent_leaving_choices_.push_back(row);
1983 equivalent_leaving_choices_.clear();
1984 current_ratio = ratio;
1985 pivot_magnitude = candidate_magnitude;
1990 if (!equivalent_leaving_choices_.empty()) {
1991 equivalent_leaving_choices_.push_back(*leaving_row);
1993 equivalent_leaving_choices_[std::uniform_int_distribution<int>(
1994 0, equivalent_leaving_choices_.size() - 1)(random_)];
2006 if (current_ratio <= 0.0) {
2010 parameters_.degenerate_ministep_factor() *
2011 parameters_.primal_feasibility_tolerance();
2012 *step_length = minimum_delta / pivot_magnitude;
2014 *step_length = current_ratio;
2021 TestPivot(entering_col, *leaving_row);
2034 if (pivot_magnitude <
2035 parameters_.small_pivot_threshold() * direction_infinity_norm_) {
2039 if (!basis_factorization_.IsRefactorized()) {
2040 VLOG(1) <<
"Refactorizing to avoid pivoting by "
2041 << direction_[*leaving_row]
2042 <<
" direction_infinity_norm_ = " << direction_infinity_norm_
2043 <<
" reduced cost = " << reduced_cost;
2044 *refactorize =
true;
2054 VLOG(1) <<
"Couldn't avoid pivoting by " << direction_[*leaving_row]
2055 <<
" direction_infinity_norm_ = " << direction_infinity_norm_
2056 <<
" reduced cost = " << reduced_cost;
2057 DCHECK_GE(std::abs(direction_[*leaving_row]),
2058 parameters_.minimum_acceptable_pivot());
2066 const bool is_reduced_cost_positive = (reduced_cost > 0.0);
2067 const bool is_leaving_coeff_positive = (direction_[*leaving_row] > 0.0);
2068 *target_bound = (is_reduced_cost_positive == is_leaving_coeff_positive)
2069 ? upper_bounds[basis_[*leaving_row]]
2070 : lower_bounds[basis_[*leaving_row]];
2075 ratio_test_stats_.leaving_choices.Add(stats_num_leaving_choices);
2076 if (!equivalent_leaving_choices_.empty()) {
2077 ratio_test_stats_.num_perfect_ties.Add(
2078 equivalent_leaving_choices_.size());
2081 ratio_test_stats_.abs_used_pivot.Add(std::abs(direction_[*leaving_row]));
2097 coeff_magnitude(_coeff_magnitude),
2098 target_bound(_target_bound) {}
2103 bool operator<(
const BreakPoint& other)
const {
2104 if (ratio == other.ratio) {
2105 if (coeff_magnitude == other.coeff_magnitude) {
2106 return row > other.row;
2108 return coeff_magnitude < other.coeff_magnitude;
2110 return ratio > other.ratio;
2121void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow(
2122 ColIndex entering_col,
Fractional reduced_cost,
bool* refactorize,
2123 RowIndex* leaving_row,
Fractional* step_length,
2130 DCHECK_NE(0.0, reduced_cost);
2131 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2132 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2136 const Fractional entering_value = variable_values_.Get(entering_col);
2137 Fractional current_ratio = (reduced_cost > 0.0)
2138 ? entering_value - lower_bounds[entering_col]
2139 : upper_bounds[entering_col] - entering_value;
2140 DCHECK_GT(current_ratio, 0.0);
2142 std::vector<BreakPoint> breakpoints;
2143 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
2144 for (
const auto e : direction_) {
2146 reduced_cost > 0.0 ? e.coefficient() : -e.coefficient();
2147 const Fractional magnitude = std::abs(direction);
2148 if (magnitude < tolerance)
continue;
2163 const ColIndex col = basis_[e.row()];
2164 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(col));
2166 const Fractional value = variable_values_.Get(col);
2167 const Fractional lower_bound = lower_bounds[col];
2168 const Fractional upper_bound = upper_bounds[col];
2169 const Fractional to_lower = (lower_bound - tolerance - value) / direction;
2170 const Fractional to_upper = (upper_bound + tolerance - value) / direction;
2174 if (to_lower >= 0.0 && to_lower < current_ratio) {
2175 breakpoints.push_back(
2176 BreakPoint(e.row(), to_lower, magnitude, lower_bound));
2178 if (to_upper >= 0.0 && to_upper < current_ratio) {
2179 breakpoints.push_back(
2180 BreakPoint(e.row(), to_upper, magnitude, upper_bound));
2186 std::make_heap(breakpoints.begin(), breakpoints.end());
2190 Fractional improvement = std::abs(reduced_cost);
2193 while (!breakpoints.empty()) {
2194 const BreakPoint top = breakpoints.front();
2202 if (top.coeff_magnitude > best_magnitude) {
2203 *leaving_row = top.row;
2204 current_ratio = top.ratio;
2205 best_magnitude = top.coeff_magnitude;
2206 *target_bound = top.target_bound;
2211 improvement -= top.coeff_magnitude;
2212 if (improvement <= 0.0)
break;
2213 std::pop_heap(breakpoints.begin(), breakpoints.end());
2214 breakpoints.pop_back();
2220 parameters_.small_pivot_threshold() * direction_infinity_norm_;
2221 if (best_magnitude < threshold && !basis_factorization_.IsRefactorized()) {
2222 *refactorize =
true;
2226 *step_length = current_ratio;
2230Status RevisedSimplex::DualChooseLeavingVariableRow(RowIndex* leaving_row,
2239 if (dual_prices_.Size() == 0) {
2240 variable_values_.RecomputeDualPrices(
2241 parameters_.dual_price_prioritize_norm());
2246 *leaving_row = dual_prices_.GetMaximum();
2249 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2250 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2251 const ColIndex leaving_col = basis_[*leaving_row];
2252 const Fractional value = variable_values_.Get(leaving_col);
2253 if (value < lower_bounds[leaving_col]) {
2254 *cost_variation = lower_bounds[leaving_col] - value;
2255 *target_bound = lower_bounds[leaving_col];
2256 DCHECK_GT(*cost_variation, 0.0);
2258 *cost_variation = upper_bounds[leaving_col] - value;
2259 *target_bound = upper_bounds[leaving_col];
2260 DCHECK_LT(*cost_variation, 0.0);
2271 if (cost == 0.0)
return false;
2281template <
bool use_dense_update>
2285 const Fractional price = dual_pricing_vector_[row];
2286 const bool is_candidate =
2287 IsDualPhaseILeavingCandidate(price, type, threshold);
2289 if (use_dense_update) {
2290 dual_prices_.DenseAddOrUpdate(row,
Square(price) / squared_norm[row]);
2292 dual_prices_.AddOrUpdate(row,
Square(price) / squared_norm[row]);
2295 dual_prices_.Remove(row);
2299void RevisedSimplex::DualPhaseIUpdatePrice(RowIndex leaving_row,
2300 ColIndex entering_col) {
2309 if (reduced_costs_.AreReducedCostsRecomputed() ||
2310 dual_edge_norms_.NeedsBasisRefactorization() ||
2311 dual_pricing_vector_.empty()) {
2316 const Fractional threshold = parameters_.ratio_test_zero_threshold();
2321 dual_edge_norms_.GetEdgeSquaredNorms();
2327 dual_pricing_vector_[leaving_row] / direction_[leaving_row];
2328 for (
const auto e : direction_) {
2329 dual_pricing_vector_[e.row()] -= e.coefficient() * step;
2330 OnDualPriceChange(squared_norms, e.row(), variable_type[basis_[e.row()]],
2333 dual_pricing_vector_[leaving_row] = step;
2337 dual_pricing_vector_[leaving_row] -=
2338 dual_infeasibility_improvement_direction_[entering_col];
2339 if (dual_infeasibility_improvement_direction_[entering_col] != 0.0) {
2340 --num_dual_infeasible_positions_;
2342 dual_infeasibility_improvement_direction_[entering_col] = 0.0;
2345 dual_infeasibility_improvement_direction_[basis_[leaving_row]] = 0.0;
2348 OnDualPriceChange(squared_norms, leaving_row, variable_type[entering_col],
2352template <
typename Cols>
2353void RevisedSimplex::DualPhaseIUpdatePriceOnReducedCostChange(
2356 bool something_to_do =
false;
2357 const DenseBitRow::ConstView can_decrease =
2358 variables_info_.GetCanDecreaseBitRow().const_view();
2359 const DenseBitRow::ConstView can_increase =
2360 variables_info_.GetCanIncreaseBitRow().const_view();
2362 const Fractional tolerance = reduced_costs_.GetDualFeasibilityTolerance();
2363 auto improvement_direction = dual_infeasibility_improvement_direction_.view();
2364 for (
const ColIndex col : cols) {
2365 const Fractional reduced_cost = reduced_costs[col];
2367 (can_increase[col] && reduced_cost < -tolerance) ? 1.0
2368 : (can_decrease[col] && reduced_cost > tolerance) ? -1.0
2370 if (sign != improvement_direction[col]) {
2372 --num_dual_infeasible_positions_;
2373 }
else if (improvement_direction[col] == 0.0) {
2374 ++num_dual_infeasible_positions_;
2376 if (!something_to_do) {
2377 initially_all_zero_scratchpad_.values.resize(num_rows_, 0.0);
2378 initially_all_zero_scratchpad_.ClearSparseMask();
2379 initially_all_zero_scratchpad_.non_zeros.clear();
2380 something_to_do =
true;
2384 num_update_price_operations_ +=
2385 10 * compact_matrix_.column(col).num_entries().value();
2386 compact_matrix_.ColumnAddMultipleToSparseScatteredColumn(
2387 col, sign - improvement_direction[col],
2388 &initially_all_zero_scratchpad_);
2389 improvement_direction[col] = sign;
2392 if (something_to_do) {
2393 initially_all_zero_scratchpad_.ClearNonZerosIfTooDense();
2394 initially_all_zero_scratchpad_.ClearSparseMask();
2396 dual_edge_norms_.GetEdgeSquaredNorms();
2399 const Fractional threshold = parameters_.ratio_test_zero_threshold();
2400 basis_factorization_.RightSolve(&initially_all_zero_scratchpad_);
2401 if (initially_all_zero_scratchpad_.non_zeros.empty()) {
2402 dual_prices_.StartDenseUpdates();
2403 for (RowIndex row(0); row < num_rows_; ++row) {
2404 if (initially_all_zero_scratchpad_[row] == 0.0)
continue;
2405 dual_pricing_vector_[row] += initially_all_zero_scratchpad_[row];
2406 OnDualPriceChange<
true>(
2407 squared_norms, row, variable_type[basis_[row]], threshold);
2409 initially_all_zero_scratchpad_.values.AssignToZero(num_rows_);
2411 for (
const auto e : initially_all_zero_scratchpad_) {
2412 dual_pricing_vector_[e.row()] += e.coefficient();
2413 OnDualPriceChange(squared_norms, e.row(),
2414 variable_type[basis_[e.row()]], threshold);
2415 initially_all_zero_scratchpad_[e.row()] = 0.0;
2418 initially_all_zero_scratchpad_.non_zeros.clear();
2422Status RevisedSimplex::DualPhaseIChooseLeavingVariableRow(
2423 RowIndex* leaving_row,
Fractional* cost_variation,
2428 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2429 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2441 if (reduced_costs_.AreReducedCostsRecomputed() ||
2442 dual_edge_norms_.NeedsBasisRefactorization() ||
2443 dual_pricing_vector_.empty()) {
2445 num_dual_infeasible_positions_ = 0;
2446 dual_pricing_vector_.AssignToZero(num_rows_);
2447 dual_prices_.ClearAndResize(num_rows_);
2448 dual_infeasibility_improvement_direction_.AssignToZero(num_cols_);
2449 DualPhaseIUpdatePriceOnReducedCostChange(
2450 variables_info_.GetIsRelevantBitRow());
2454 DualPhaseIUpdatePriceOnReducedCostChange(update_row_.GetNonZeroPositions());
2459 if (num_dual_infeasible_positions_ == 0)
return Status::OK();
2461 *leaving_row = dual_prices_.GetMaximum();
2466 *cost_variation = dual_pricing_vector_[*leaving_row];
2467 const ColIndex leaving_col = basis_[*leaving_row];
2468 if (*cost_variation < 0.0) {
2469 *target_bound = upper_bounds[leaving_col];
2471 *target_bound = lower_bounds[leaving_col];
2477template <
typename BoxedVariableCols>
2478void RevisedSimplex::MakeBoxedVariableDualFeasible(
2479 const BoxedVariableCols& cols,
bool update_basic_values) {
2481 std::vector<ColIndex> changed_cols;
2491 const Fractional threshold = reduced_costs_.GetDualFeasibilityTolerance();
2493 const DenseRow& variable_values = variable_values_.GetDenseRow();
2495 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2496 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2498 for (
const ColIndex col : cols) {
2499 const Fractional reduced_cost = reduced_costs[col];
2501 DCHECK(variables_info_.GetTypeRow()[col] ==
2504 DCHECK(variable_values[col] == lower_bounds[col] ||
2505 variable_values[col] == upper_bounds[col] ||
2508 variables_info_.UpdateToNonBasicStatus(col,
2510 changed_cols.push_back(col);
2511 }
else if (reduced_cost < -threshold &&
2513 variables_info_.UpdateToNonBasicStatus(col,
2515 changed_cols.push_back(col);
2519 if (!changed_cols.empty()) {
2520 iteration_stats_.num_dual_flips.Add(changed_cols.size());
2521 variable_values_.UpdateGivenNonBasicVariables(changed_cols,
2522 update_basic_values);
2526Fractional RevisedSimplex::ComputeStepToMoveBasicVariableToBound(
2527 RowIndex leaving_row,
Fractional target_bound) {
2531 const ColIndex leaving_col = basis_[leaving_row];
2532 const Fractional leaving_variable_value = variable_values_.Get(leaving_col);
2533 Fractional unscaled_step = leaving_variable_value - target_bound;
2542 return unscaled_step / direction_[leaving_row];
2545bool RevisedSimplex::TestPivot(ColIndex entering_col, RowIndex leaving_row) {
2546 VLOG(1) <<
"Test pivot.";
2548 const ColIndex leaving_col = basis_[leaving_row];
2549 basis_[leaving_row] = entering_col;
2553 CompactSparseMatrixView basis_matrix(&compact_matrix_, &basis_);
2554 const bool is_ok = test_lu_.ComputeFactorization(basis_matrix).ok();
2555 basis_[leaving_row] = leaving_col;
2562void RevisedSimplex::PermuteBasis() {
2568 basis_factorization_.GetColumnPermutation();
2569 if (col_perm.empty())
return;
2576 if (!dual_pricing_vector_.empty()) {
2580 &dual_pricing_vector_,
2581 &tmp_dual_pricing_vector_);
2585 reduced_costs_.UpdateDataOnBasisPermutation();
2586 dual_edge_norms_.UpdateDataOnBasisPermutation(col_perm);
2590 basis_factorization_.SetColumnPermutationToIdentity();
2593Status RevisedSimplex::UpdateAndPivot(ColIndex entering_col,
2594 RowIndex leaving_row,
2610 if (update_row_.IsComputedFor(leaving_row)) {
2611 pivot_from_update_row = update_row_.GetCoefficient(entering_col);
2615 update_row_.ComputeUnitRowLeftInverse(leaving_row);
2616 pivot_from_update_row = compact_matrix_.ColumnScalarProduct(
2617 entering_col, update_row_.GetUnitRowLeftInverse().values);
2620 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2621 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2622 const ColIndex leaving_col = basis_[leaving_row];
2624 lower_bounds[leaving_col] == upper_bounds[leaving_col]
2626 : target_bound == lower_bounds[leaving_col]
2629 if (variable_values_.Get(leaving_col) != target_bound) {
2630 ratio_test_stats_.bound_shift.Add(variable_values_.Get(leaving_col) -
2633 UpdateBasis(entering_col, leaving_row, leaving_variable_status);
2636 const Fractional pivot_from_direction = direction_[leaving_row];
2638 std::abs(pivot_from_update_row - pivot_from_direction);
2639 if (diff > parameters_.refactorization_threshold() *
2640 (1.0 + std::min(std::abs(pivot_from_update_row),
2641 std::abs(pivot_from_direction)))) {
2642 VLOG(1) <<
"Refactorizing: imprecise pivot " << pivot_from_direction
2643 <<
" diff = " << diff;
2646 if (basis_factorization_.NumUpdates() < 10) {
2647 Fractional threshold = parameters_.lu_factorization_pivot_threshold();
2648 threshold = std::min(threshold * 1.5, 0.9);
2649 VLOG(1) <<
"Increasing LU pivot threshold " << threshold;
2650 parameters_.set_lu_factorization_pivot_threshold(threshold);
2651 basis_factorization_.SetParameters(parameters_);
2654 last_refactorization_reason_ = RefactorizationReason::IMPRECISE_PIVOT;
2658 basis_factorization_.Update(entering_col, leaving_row, direction_));
2660 if (basis_factorization_.IsRefactorized()) {
2666Status RevisedSimplex::RefactorizeBasisIfNeeded(
bool* refactorize) {
2668 if (*refactorize && !basis_factorization_.IsRefactorized()) {
2670 update_row_.Invalidate();
2673 *refactorize =
false;
2678 if (col >= integrality_scale_.size()) {
2679 integrality_scale_.resize(col + 1, 0.0);
2681 integrality_scale_[col] = scale;
2686 Cleanup update_deterministic_time_on_return(
2693 std::vector<ColIndex> candidates;
2696 if (std::abs(rc[col]) < 1e-9) candidates.push_back(col);
2699 bool refactorize =
false;
2702 for (
int i = 0; i < 10; ++i) {
2705 if (num_pivots >= 5)
break;
2706 if (candidates.empty())
break;
2710 std::uniform_int_distribution<int>(0, candidates.size() - 1)(random_);
2711 const ColIndex entering_col = candidates[index];
2712 std::swap(candidates[index], candidates.back());
2713 candidates.pop_back();
2723 if (reduced_costs_.NeedsBasisRefactorization()) refactorize =
true;
2727 ComputeDirection(entering_col);
2729 RowIndex leaving_row;
2731 bool local_refactorize =
false;
2733 ChooseLeavingVariableRow(entering_col, fake_rc, &local_refactorize,
2734 &leaving_row, &step_length, &target_bound));
2736 if (local_refactorize)
continue;
2738 if (std::abs(step_length) <= 1e-6)
continue;
2739 if (leaving_row !=
kInvalidRow && std::abs(direction_[leaving_row]) < 0.1) {
2742 const Fractional step = (fake_rc > 0.0) ? -step_length : step_length;
2748 const auto get_diff = [
this](ColIndex col,
Fractional old_value,
2750 if (col >= integrality_scale_.size() || integrality_scale_[col] == 0.0) {
2753 const Fractional s = integrality_scale_[col];
2754 return (std::abs(new_value * s - std::round(new_value * s)) -
2755 std::abs(old_value * s - std::round(old_value * s)));
2757 Fractional diff = get_diff(entering_col, variable_values_.Get(entering_col),
2758 variable_values_.Get(entering_col) + step);
2759 for (
const auto e : direction_) {
2760 const ColIndex col = basis_[e.row()];
2761 const Fractional old_value = variable_values_.Get(col);
2762 const Fractional new_value = old_value - e.coefficient() * step;
2763 diff += get_diff(col, old_value, new_value);
2767 if (diff > -1e-2)
continue;
2772 variable_values_.UpdateOnPivoting(direction_, entering_col, step);
2777 SetNonBasicVariableStatusAndDeriveValue(entering_col,
2779 }
else if (step < 0.0) {
2780 SetNonBasicVariableStatusAndDeriveValue(entering_col,
2787 const ColIndex leaving_col = basis_[leaving_row];
2788 update_row_.ComputeUpdateRow(leaving_row);
2797 primal_edge_norms_.UpdateBeforeBasisPivot(
2798 entering_col, leaving_col, leaving_row, direction_, &update_row_);
2799 dual_edge_norms_.UpdateBeforeBasisPivot(
2800 entering_col, leaving_row, direction_,
2801 update_row_.GetUnitRowLeftInverse());
2806 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row, direction_,
2809 const Fractional dir = -direction_[leaving_row] * step;
2810 const bool is_degenerate =
2812 (dir > 0.0 && variable_values_.Get(leaving_col) >= target_bound) ||
2813 (dir < 0.0 && variable_values_.Get(leaving_col) <= target_bound);
2814 if (!is_degenerate) {
2815 variable_values_.Set(leaving_col, target_bound);
2818 UpdateAndPivot(entering_col, leaving_row, target_bound));
2821 VLOG(1) <<
"Polish num_pivots: " << num_pivots <<
" gain:" << total_gain;
2842 Cleanup update_deterministic_time_on_return(
2844 num_consecutive_degenerate_iterations_ = 0;
2845 bool refactorize =
false;
2846 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
2850 primal_prices_.ForceRecomputation();
2852 if (phase_ == Phase::FEASIBILITY) {
2855 objective_.AssignToZero(num_cols_);
2856 variable_values_.UpdatePrimalPhaseICosts(
2857 util::IntegerRange<RowIndex>(RowIndex(0), num_rows_), &objective_);
2858 reduced_costs_.ResetForNewObjective();
2870 if (!refactorize && reduced_costs_.NeedsBasisRefactorization()) {
2871 last_refactorization_reason_ = RefactorizationReason::RC;
2874 if (!refactorize && primal_edge_norms_.NeedsBasisRefactorization()) {
2875 last_refactorization_reason_ = RefactorizationReason::NORM;
2880 if (basis_factorization_.IsRefactorized()) {
2881 CorrectErrorsOnVariableValues();
2882 DisplayIterationInfo(
true, last_refactorization_reason_);
2883 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
2885 if (phase_ == Phase::FEASIBILITY) {
2888 if (variable_values_.UpdatePrimalPhaseICosts(
2889 util::IntegerRange<RowIndex>(RowIndex(0), num_rows_),
2891 reduced_costs_.ResetForNewObjective();
2897 if (phase_ == Phase::OPTIMIZATION &&
2898 ComputeObjectiveValue() < primal_objective_limit_) {
2899 VLOG(1) <<
"Stopping the primal simplex because"
2900 <<
" the objective limit " << primal_objective_limit_
2901 <<
" has been reached.";
2903 objective_limit_reached_ =
true;
2906 }
else if (phase_ == Phase::FEASIBILITY) {
2909 if (variable_values_.UpdatePrimalPhaseICosts(direction_.non_zeros,
2911 reduced_costs_.ResetForNewObjective();
2915 const ColIndex entering_col = primal_prices_.GetBestEnteringColumn();
2917 if (reduced_costs_.AreReducedCostsPrecise() &&
2918 basis_factorization_.IsRefactorized()) {
2919 if (phase_ == Phase::FEASIBILITY) {
2921 variable_values_.ComputeMaximumPrimalInfeasibility();
2922 if (primal_infeasibility <
2923 parameters_.primal_feasibility_tolerance()) {
2926 VLOG(1) <<
"Infeasible problem! infeasibility = "
2927 << primal_infeasibility;
2936 VLOG(1) <<
"Optimal reached, double checking...";
2937 reduced_costs_.MakeReducedCostsPrecise();
2939 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
2943 DCHECK(reduced_costs_.IsValidPrimalEnteringCandidate(entering_col));
2946 ComputeDirection(entering_col);
2951 if (!primal_edge_norms_.TestEnteringEdgeNormPrecision(entering_col,
2953 primal_prices_.RecomputePriceAt(entering_col);
2957 reduced_costs_.TestEnteringReducedCostPrecision(entering_col,
2962 primal_prices_.RecomputePriceAt(entering_col);
2963 if (!reduced_costs_.IsValidPrimalEnteringCandidate(entering_col)) {
2964 reduced_costs_.MakeReducedCostsPrecise();
2965 VLOG(1) <<
"Skipping col #" << entering_col
2966 <<
" whose reduced cost is no longer valid under precise reduced "
2976 if (num_iterations_ == parameters_.max_number_of_iterations())
break;
2979 RowIndex leaving_row;
2981 if (phase_ == Phase::FEASIBILITY) {
2982 PrimalPhaseIChooseLeavingVariableRow(entering_col, reduced_cost,
2983 &refactorize, &leaving_row,
2984 &step_length, &target_bound);
2987 ChooseLeavingVariableRow(entering_col, reduced_cost, &refactorize,
2988 &leaving_row, &step_length, &target_bound));
2991 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
2999 if (!basis_factorization_.IsRefactorized() ||
3000 !reduced_costs_.AreReducedCostsPrecise()) {
3001 VLOG(1) <<
"Infinite step length, double checking...";
3002 reduced_costs_.MakeReducedCostsPrecise();
3004 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3007 if (phase_ == Phase::FEASIBILITY) {
3009 VLOG(1) <<
"Unbounded feasibility problem !?";
3013 solution_primal_ray_.AssignToZero(num_cols_);
3014 for (RowIndex row(0); row < num_rows_; ++row) {
3015 const ColIndex col = basis_[row];
3016 solution_primal_ray_[col] = -direction_[row];
3018 solution_primal_ray_[entering_col] = 1.0;
3019 if (reduced_cost > 0.0) {
3026 Fractional step = (reduced_cost > 0.0) ? -step_length : step_length;
3027 if (phase_ == Phase::FEASIBILITY && leaving_row !=
kInvalidRow) {
3037 step = ComputeStepToMoveBasicVariableToBound(leaving_row, target_bound);
3041 const ColIndex leaving_col =
3047 bool is_degenerate =
false;
3049 Fractional dir = -direction_[leaving_row] * step;
3052 (dir > 0.0 && variable_values_.Get(leaving_col) >= target_bound) ||
3053 (dir < 0.0 && variable_values_.Get(leaving_col) <= target_bound);
3057 if (!is_degenerate) {
3058 DCHECK_EQ(step, ComputeStepToMoveBasicVariableToBound(leaving_row,
3063 variable_values_.UpdateOnPivoting(direction_, entering_col, step);
3066 primal_edge_norms_.UpdateBeforeBasisPivot(
3067 entering_col, basis_[leaving_row], leaving_row, direction_,
3069 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row,
3070 direction_, &update_row_);
3071 primal_prices_.UpdateBeforeBasisPivot(entering_col, &update_row_);
3072 if (!is_degenerate) {
3078 variable_values_.Set(leaving_col, target_bound);
3081 UpdateAndPivot(entering_col, leaving_row, target_bound));
3083 if (is_degenerate) {
3084 timer.AlsoUpdate(&iteration_stats_.degenerate);
3086 timer.AlsoUpdate(&iteration_stats_.normal);
3093 variables_info_.GetTypeRow()[entering_col]);
3095 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3097 }
else if (step < 0.0) {
3098 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3101 primal_prices_.SetAndDebugCheckThatColumnIsDualFeasible(entering_col);
3105 if (phase_ == Phase::FEASIBILITY && leaving_row !=
kInvalidRow) {
3107 variable_values_.SetNonBasicVariableValueFromStatus(leaving_col);
3110 reduced_costs_.SetNonBasicVariableCostToZero(leaving_col,
3111 &objective_[leaving_col]);
3112 primal_prices_.RecomputePriceAt(leaving_col);
3116 if (step_length == 0.0) {
3117 num_consecutive_degenerate_iterations_++;
3119 if (num_consecutive_degenerate_iterations_ > 0) {
3120 iteration_stats_.degenerate_run_size.Add(
3121 num_consecutive_degenerate_iterations_);
3122 num_consecutive_degenerate_iterations_ = 0;
3127 if (num_consecutive_degenerate_iterations_ > 0) {
3128 iteration_stats_.degenerate_run_size.Add(
3129 num_consecutive_degenerate_iterations_);
3145Status RevisedSimplex::DualMinimize(
bool feasibility_phase,
3147 Cleanup update_deterministic_time_on_return(
3149 num_consecutive_degenerate_iterations_ = 0;
3150 bool refactorize =
false;
3151 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
3153 bound_flip_candidates_.clear();
3156 RowIndex leaving_row;
3161 ColIndex entering_col;
3174 const bool old_refactorize_value = refactorize;
3175 if (!refactorize && reduced_costs_.NeedsBasisRefactorization()) {
3176 last_refactorization_reason_ = RefactorizationReason::RC;
3179 if (!refactorize && dual_edge_norms_.NeedsBasisRefactorization()) {
3180 last_refactorization_reason_ = RefactorizationReason::NORM;
3187 if (basis_factorization_.IsRefactorized()) {
3201 if (feasibility_phase || old_refactorize_value) {
3202 reduced_costs_.MakeReducedCostsPrecise();
3213 if (!feasibility_phase) {
3214 MakeBoxedVariableDualFeasible(
3215 variables_info_.GetNonBasicBoxedVariables(),
3217 variable_values_.RecomputeBasicVariableValues();
3218 variable_values_.RecomputeDualPrices(
3219 parameters_.dual_price_prioritize_norm());
3228 if (phase_ == Phase::OPTIMIZATION &&
3230 ComputeObjectiveValue() > dual_objective_limit_) {
3232 "Stopping the dual simplex because"
3233 " the objective limit ",
3234 dual_objective_limit_,
" has been reached.");
3236 objective_limit_reached_ =
true;
3241 DisplayIterationInfo(
false, last_refactorization_reason_);
3242 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
3246 if (!feasibility_phase) {
3249 MakeBoxedVariableDualFeasible(bound_flip_candidates_,
3251 bound_flip_candidates_.clear();
3255 variable_values_.UpdateDualPrices(direction_.non_zeros);
3259 if (feasibility_phase) {
3261 &leaving_row, &cost_variation, &target_bound));
3264 &leaving_row, &cost_variation, &target_bound));
3269 if (!basis_factorization_.IsRefactorized() ||
3270 reduced_costs_.HasCostShift()) {
3271 VLOG(1) <<
"Optimal reached, double checking.";
3272 reduced_costs_.ClearAndRemoveCostShifts();
3275 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3278 if (feasibility_phase) {
3283 if (num_dual_infeasible_positions_ == 0) {
3286 VLOG(1) <<
"DUAL infeasible in dual phase I.";
3296 update_row_.ComputeUnitRowLeftInverse(leaving_row);
3297 if (!dual_edge_norms_.TestPrecision(leaving_row,
3298 update_row_.GetUnitRowLeftInverse())) {
3302 if (feasibility_phase) {
3303 const Fractional price = dual_pricing_vector_[leaving_row];
3305 dual_edge_norms_.GetEdgeSquaredNorms();
3306 dual_prices_.AddOrUpdate(leaving_row,
3307 Square(price) / squared_norms[leaving_row]);
3309 variable_values_.UpdateDualPrices({leaving_row});
3313 update_row_.ComputeUpdateRow(leaving_row);
3315 if (feasibility_phase) {
3317 reduced_costs_.AreReducedCostsPrecise(), update_row_, cost_variation,
3321 reduced_costs_.AreReducedCostsPrecise(), update_row_, cost_variation,
3322 &bound_flip_candidates_, &entering_col));
3327 if (!reduced_costs_.AreReducedCostsPrecise()) {
3328 VLOG(1) <<
"No entering column. Double checking...";
3331 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3334 DCHECK(basis_factorization_.IsRefactorized());
3335 if (feasibility_phase) {
3337 VLOG(1) <<
"Unbounded dual feasibility problem !?";
3341 solution_dual_ray_ =
3342 Transpose(update_row_.GetUnitRowLeftInverse().values);
3343 update_row_.ComputeFullUpdateRow(leaving_row,
3344 &solution_dual_ray_row_combination_);
3345 if (cost_variation < 0) {
3347 ChangeSign(&solution_dual_ray_row_combination_);
3357 const Fractional entering_coeff = update_row_.GetCoefficient(entering_col);
3358 if (std::abs(entering_coeff) < parameters_.dual_small_pivot_threshold() &&
3359 !reduced_costs_.AreReducedCostsPrecise()) {
3360 VLOG(1) <<
"Trying not to pivot by " << entering_coeff;
3363 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
3367 ComputeDirection(entering_col);
3373 if (std::abs(direction_[leaving_row]) <
3374 parameters_.small_pivot_threshold() * direction_infinity_norm_) {
3375 if (!reduced_costs_.AreReducedCostsPrecise()) {
3376 VLOG(1) <<
"Trying not pivot by " << entering_coeff <<
" ("
3377 << direction_[leaving_row]
3378 <<
") because the direction has a norm of "
3379 << direction_infinity_norm_;
3382 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
3391 if (std::abs(direction_[leaving_row]) <= 1e-20) {
3392 const std::string error_message = absl::StrCat(
3393 "trying to pivot with number too small: ", direction_[leaving_row]);
3402 if (num_iterations_ == parameters_.max_number_of_iterations()) {
3416 const bool increasing_rc_is_needed =
3417 (cost_variation > 0.0) == (entering_coeff > 0.0);
3418 reduced_costs_.ShiftCostIfNeeded(increasing_rc_is_needed, entering_col);
3421 if (reduced_costs_.StepIsDualDegenerate(increasing_rc_is_needed,
3423 timer.AlsoUpdate(&iteration_stats_.degenerate);
3425 timer.AlsoUpdate(&iteration_stats_.normal);
3434 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row, direction_,
3436 dual_edge_norms_.UpdateBeforeBasisPivot(
3437 entering_col, leaving_row, direction_,
3438 update_row_.GetUnitRowLeftInverse());
3443 if (feasibility_phase) {
3444 DualPhaseIUpdatePrice(leaving_row, entering_col);
3447 ComputeStepToMoveBasicVariableToBound(leaving_row, target_bound);
3448 variable_values_.UpdateOnPivoting(direction_, entering_col, primal_step);
3452 const ColIndex leaving_col = basis_[leaving_row];
3454 UpdateAndPivot(entering_col, leaving_row, target_bound));
3460 variable_values_.SetNonBasicVariableValueFromStatus(leaving_col);
3469 Cleanup update_deterministic_time_on_return(
3471 bool refactorize =
false;
3475 primal_edge_norms_.Clear();
3476 dual_edge_norms_.Clear();
3477 update_row_.Invalidate();
3478 reduced_costs_.ClearAndRemoveCostShifts();
3480 std::vector<ColIndex> super_basic_cols;
3481 for (
const ColIndex col : variables_info_.GetNotBasicBitRow()) {
3483 variable_values_.Get(col) != 0) {
3484 super_basic_cols.push_back(col);
3488 while (!super_basic_cols.empty()) {
3494 if (basis_factorization_.IsRefactorized()) {
3495 CorrectErrorsOnVariableValues();
3496 DisplayIterationInfo(
true);
3500 ColIndex entering_col = super_basic_cols.back();
3502 DCHECK(variables_info_.GetCanDecreaseBitRow()[entering_col]);
3503 DCHECK(variables_info_.GetCanIncreaseBitRow()[entering_col]);
3514 const Fractional entering_value = variable_values_.Get(entering_col);
3515 if (variables_info_.GetTypeRow()[entering_col] ==
3517 if (entering_value > 0) {
3524 variables_info_.GetVariableUpperBounds()[entering_col] -
3528 variables_info_.GetVariableLowerBounds()[entering_col];
3529 if (diff_lb <= diff_ub) {
3537 ComputeDirection(entering_col);
3540 RowIndex leaving_row;
3544 &refactorize, &leaving_row,
3545 &step_length, &target_bound));
3547 if (refactorize)
continue;
3550 super_basic_cols.pop_back();
3553 if (variables_info_.GetTypeRow()[entering_col] ==
3555 step_length = std::fabs(entering_value);
3557 VLOG(1) <<
"Infinite step for bounded variable ?!";
3563 const Fractional step = (fake_rc > 0.0) ? -step_length : step_length;
3566 const ColIndex leaving_col =
3575 bool is_degenerate =
false;
3577 Fractional dir = -direction_[leaving_row] * step;
3580 (dir > 0.0 && variable_values_.Get(leaving_col) >= target_bound) ||
3581 (dir < 0.0 && variable_values_.Get(leaving_col) <= target_bound);
3585 if (!is_degenerate) {
3586 DCHECK_EQ(step, ComputeStepToMoveBasicVariableToBound(leaving_row,
3591 variable_values_.UpdateOnPivoting(direction_, entering_col, step);
3593 if (!is_degenerate) {
3599 variable_values_.Set(leaving_col, target_bound);
3602 UpdateAndPivot(entering_col, leaving_row, target_bound));
3604 if (is_degenerate) {
3605 timer.AlsoUpdate(&iteration_stats_.degenerate);
3607 timer.AlsoUpdate(&iteration_stats_.normal);
3614 if (variables_info_.GetTypeRow()[entering_col] ==
3616 variable_values_.Set(entering_col, 0.0);
3617 }
else if (step > 0.0) {
3618 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3620 }
else if (step < 0.0) {
3621 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3630 if (!super_basic_cols.empty()) {
3631 SOLVER_LOG(logger_,
"Push terminated early with ", super_basic_cols.size(),
3632 " super-basic variables remaining.");
3642ColIndex RevisedSimplex::SlackColIndex(RowIndex row)
const {
3649 result.append(iteration_stats_.StatString());
3650 result.append(ratio_test_stats_.StatString());
3651 result.append(entering_variable_.StatString());
3652 result.append(dual_prices_.StatString());
3653 result.append(reduced_costs_.StatString());
3654 result.append(variable_values_.StatString());
3655 result.append(primal_edge_norms_.StatString());
3656 result.append(dual_edge_norms_.StatString());
3657 result.append(update_row_.StatString());
3658 result.append(basis_factorization_.StatString());
3659 result.append(function_stats_.StatString());
3663void RevisedSimplex::DisplayAllStats() {
3664 if (absl::GetFlag(FLAGS_simplex_display_stats)) {
3666 absl::FPrintF(stderr,
"%s", GetPrettySolverStats());
3670Fractional RevisedSimplex::ComputeObjectiveValue()
const {
3676Fractional RevisedSimplex::ComputeInitialProblemObjectiveValue()
const {
3679 objective_,
Transpose(variable_values_.GetDenseRow()));
3680 return objective_scaling_factor_ * (sum + objective_offset_);
3685 dual_prices_.SetRandom(random_);
3686 entering_variable_.SetRandom(random_);
3687 reduced_costs_.SetRandom(random_);
3688 primal_prices_.SetRandom(random_);
3693 deterministic_random_.seed(parameters.
random_seed());
3694 absl_random_ = absl::BitGen(absl::SeedSeq({parameters.
random_seed()}));
3696 ? absl::BitGenRef(absl_random_)
3697 : absl::BitGenRef(deterministic_random_)));
3698 initial_parameters_ = parameters;
3699 parameters_ = parameters;
3700 PropagateParameters();
3703void RevisedSimplex::PropagateParameters() {
3713void RevisedSimplex::DisplayIterationInfo(
bool primal,
3714 RefactorizationReason reason) {
3716 const std::string first_word = primal ?
"Primal " :
"Dual ";
3723 if (reason != RefactorizationReason::DEFAULT) {
3725 case RefactorizationReason::DEFAULT:
3726 info =
" [default]";
3728 case RefactorizationReason::SMALL_PIVOT:
3729 info =
" [small pivot]";
3731 case RefactorizationReason::IMPRECISE_PIVOT:
3732 info =
" [imprecise pivot]";
3734 case RefactorizationReason::NORM:
3737 case RefactorizationReason::RC:
3738 info =
" [reduced costs]";
3740 case RefactorizationReason::VAR_VALUES:
3741 info =
" [var values]";
3743 case RefactorizationReason::FINAL_CHECK:
3750 case Phase::FEASIBILITY: {
3751 const int64_t iter = num_iterations_;
3754 if (parameters_.use_dual_simplex()) {
3755 if (parameters_.use_dedicated_dual_feasibility_algorithm()) {
3756 objective = reduced_costs_.ComputeSumOfDualInfeasibilities();
3761 objective_,
Transpose(variable_values_.GetDenseRow()));
3763 name =
"sum_dual_infeasibilities";
3765 objective = variable_values_.ComputeSumOfPrimalInfeasibilities();
3766 name =
"sum_primal_infeasibilities";
3769 SOLVER_LOG(logger_, first_word,
"feasibility phase, iteration # ", iter,
3770 ", ", name,
" = ", absl::StrFormat(
"%.15E", objective), info);
3773 case Phase::OPTIMIZATION: {
3774 const int64_t iter = num_iterations_ - num_feasibility_iterations_;
3780 const Fractional objective = ComputeInitialProblemObjectiveValue();
3781 SOLVER_LOG(logger_, first_word,
"optimization phase, iteration # ", iter,
3782 ", objective = ", absl::StrFormat(
"%.15E", objective), info);
3786 const int64_t iter = num_iterations_ - num_feasibility_iterations_ -
3787 num_optimization_iterations_;
3788 SOLVER_LOG(logger_, first_word,
"push phase, iteration # ", iter,
3789 ", remaining_variables_to_push = ",
3790 ComputeNumberOfSuperBasicVariables(), info);
3795void RevisedSimplex::DisplayErrors() {
3796 if (!logger_->LoggingIsEnabled())
return;
3799 SOLVER_LOG(logger_,
"Primal infeasibility (bounds) = ",
3800 variable_values_.ComputeMaximumPrimalInfeasibility());
3801 SOLVER_LOG(logger_,
"Primal residual |A.x - b| = ",
3802 variable_values_.ComputeMaximumPrimalResidual());
3803 SOLVER_LOG(logger_,
"Dual infeasibility (reduced costs) = ",
3804 reduced_costs_.ComputeMaximumDualInfeasibility());
3805 SOLVER_LOG(logger_,
"Dual residual |c_B - y.B| = ",
3806 reduced_costs_.ComputeMaximumDualResidual());
3811std::string StringifyMonomialWithFlags(
const Fractional a,
3812 absl::string_view x) {
3814 a, x, absl::GetFlag(FLAGS_simplex_display_numbers_as_fractions));
3820std::string StringifyWithFlags(
const Fractional x) {
3822 absl::GetFlag(FLAGS_simplex_display_numbers_as_fractions));
3827std::string RevisedSimplex::SimpleVariableInfo(ColIndex col)
const {
3829 VariableType variable_type = variables_info_.GetTypeRow()[col];
3830 VariableStatus variable_status = variables_info_.GetStatusRow()[col];
3831 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
3832 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
3833 absl::StrAppendFormat(&output,
"%d (%s) = %s, %s, %s, [%s,%s]", col.value(),
3834 variable_name_[col],
3835 StringifyWithFlags(variable_values_.Get(col)),
3838 StringifyWithFlags(lower_bounds[col]),
3839 StringifyWithFlags(upper_bounds[col]));
3843void RevisedSimplex::DisplayInfoOnVariables()
const {
3844 if (VLOG_IS_ON(3)) {
3845 for (ColIndex col(0); col < num_cols_; ++col) {
3846 const Fractional variable_value = variable_values_.Get(col);
3847 const Fractional objective_coefficient = objective_[col];
3849 objective_coefficient * variable_value;
3850 VLOG(3) << SimpleVariableInfo(col) <<
". " << variable_name_[col] <<
" = "
3851 << StringifyWithFlags(variable_value) <<
" * "
3852 << StringifyWithFlags(objective_coefficient)
3853 <<
"(obj) = " << StringifyWithFlags(objective_contribution);
3855 VLOG(3) <<
"------";
3859void RevisedSimplex::DisplayVariableBounds() {
3860 if (VLOG_IS_ON(3)) {
3862 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
3863 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
3864 for (ColIndex col(0); col < num_cols_; ++col) {
3865 switch (variable_type[col]) {
3869 VLOG(3) << variable_name_[col]
3870 <<
" >= " << StringifyWithFlags(lower_bounds[col]) <<
";";
3873 VLOG(3) << variable_name_[col]
3874 <<
" <= " << StringifyWithFlags(upper_bounds[col]) <<
";";
3877 VLOG(3) << StringifyWithFlags(lower_bounds[col])
3878 <<
" <= " << variable_name_[col]
3879 <<
" <= " << StringifyWithFlags(upper_bounds[col]) <<
";";
3882 VLOG(3) << variable_name_[col] <<
" = "
3883 << StringifyWithFlags(lower_bounds[col]) <<
";";
3886 LOG(DFATAL) <<
"Column " << col <<
" has no meaningful status.";
3893util_intops::StrongVector<RowIndex, SparseRow>
3896 for (ColIndex col(0); col < num_cols_; ++col) {
3897 ComputeDirection(col);
3898 for (
const auto e : direction_) {
3899 if (column_scales ==
nullptr) {
3900 dictionary[e.row()].SetCoefficient(col, e.coefficient());
3904 col < column_scales->
size() ? (*column_scales)[col] : 1.0;
3906 ? (*column_scales)[
GetBasis(e.row())]
3908 dictionary[e.row()].SetCoefficient(
3909 col, direction_[e.row()] * (numerator / denominator));
3918 Status status = Initialize(linear_program);
3920 variable_values_.RecomputeBasicVariableValues();
3921 solution_objective_value_ = ComputeInitialProblemObjectiveValue();
3925void RevisedSimplex::DisplayRevisedSimplexDebugInfo() {
3926 if (VLOG_IS_ON(3)) {
3927 variable_name_.
resize(num_cols_,
"");
3928 DisplayInfoOnVariables();
3930 std::string output =
"z = " + StringifyWithFlags(ComputeObjectiveValue());
3933 absl::StrAppend(&output, StringifyMonomialWithFlags(reduced_costs[col],
3934 variable_name_[col]));
3936 VLOG(3) << output <<
";";
3938 const RevisedSimplexDictionary dictionary(
nullptr,
this);
3940 for (
const SparseRow& row : dictionary) {
3942 ColIndex basic_col = basis_[r];
3943 absl::StrAppend(&output, variable_name_[basic_col],
" = ",
3944 StringifyWithFlags(variable_values_.
Get(basic_col)));
3945 for (
const SparseRowEntry e : row) {
3946 if (e.col() != basic_col) {
3947 absl::StrAppend(&output,
3948 StringifyMonomialWithFlags(e.coefficient(),
3949 variable_name_[e.col()]));
3952 VLOG(3) << output <<
";";
3954 VLOG(3) <<
"------";
3955 DisplayVariableBounds();
3960void RevisedSimplex::DisplayProblem() {
3963 if (VLOG_IS_ON(3)) {
3964 variable_name_.resize(num_cols_,
"");
3965 DisplayInfoOnVariables();
3966 std::string output =
"min: ";
3967 bool has_objective =
false;
3968 for (ColIndex col(0); col < num_cols_; ++col) {
3970 has_objective |= (coeff != 0.0);
3971 absl::StrAppend(&output,
3972 StringifyMonomialWithFlags(coeff, variable_name_[col]));
3974 if (!has_objective) {
3975 absl::StrAppend(&output,
" 0");
3977 VLOG(3) << output <<
";";
3978 for (RowIndex row(0); row < num_rows_; ++row) {
3980 for (ColIndex col(0); col < num_cols_; ++col) {
3981 absl::StrAppend(&output,
3982 StringifyMonomialWithFlags(
3983 compact_matrix_.column(col).LookUpCoefficient(row),
3984 variable_name_[col]));
3986 VLOG(3) << output <<
" = 0;";
3988 VLOG(3) <<
"------";
3992void RevisedSimplex::AdvanceDeterministicTime(TimeLimit*
time_limit) {
3995 const double deterministic_time_delta =
3996 current_deterministic_time - last_deterministic_time_update_;
3997 time_limit->AdvanceDeterministicTime(deterministic_time_delta);
3998 last_deterministic_time_update_ = current_deterministic_time;
4001#undef DCHECK_COL_BOUNDS
4002#undef DCHECK_ROW_BOUNDS
bool LoggingIsEnabled() const
Returns true iff logging is enabled.
void SetParameters(const GlopParameters ¶meters)
Sets the parameters for this component.
void SetTimeLimit(TimeLimit *time_limit)
void SetParameters(const GlopParameters ¶meters)
Sets the algorithm parameters.
void SetParameters(const GlopParameters ¶meters)
Sets the parameters.
bool use_absl_random() const
static constexpr InitialBasisHeuristic NONE
static constexpr InitialBasisHeuristic BIXBY
::int32_t random_seed() const
static constexpr InitialBasisHeuristic TRIANGULAR
static constexpr InitialBasisHeuristic MAROS
const DenseRow & objective_coefficients() const
Returns the objective coefficients (or cost) of variables as a row vector.
bool IsMaximizationProblem() const
void SetTimeLimit(TimeLimit *time_limit)
void SetParameters(const GlopParameters ¶meters)
Sets the algorithm parameters.
void SetParameters(const GlopParameters ¶meters)
Sets the pricing parameters. This does not change the pricing rule.
DenseRow::ConstView GetReducedCosts()
Fractional GetVariableValue(ColIndex col) const
const DenseColumn & GetDualRay() const
void SetIntegralityScale(ColIndex col, Fractional scale)
const DenseRow & GetReducedCosts() const
ColIndex GetBasis(RowIndex row) const
ProblemStatus GetProblemStatus() const
ConstraintStatus GetConstraintStatus(RowIndex row) const
Fractional GetObjectiveValue() const
const BasisState & GetState() const
void LoadStateForNextSolve(const BasisState &state)
Uses the given state as a warm-start for the next Solve() call.
Fractional GetDualValue(RowIndex row) const
Fractional GetReducedCost(ColIndex col) const
ColIndex GetProblemNumCols() const
void ClearStateForNextSolve()
Fractional GetConstraintActivity(RowIndex row) const
RowMajorSparseMatrix ComputeDictionary(const DenseRow *column_scales)
const DenseRow & GetDualRayRowCombination() const
This is the "dual ray" linear combination of the matrix rows.
RowIndex GetProblemNumRows() const
Getters to retrieve all the information computed by the last Solve().
void ComputeBasicVariablesForState(const LinearProgram &linear_program, const BasisState &state)
const DenseRow & GetPrimalRay() const
std::string StatString()
Returns statistics about this class as a string.
const BasisFactorization & GetBasisFactorization() const
VariableStatus GetVariableStatus(ColIndex col) const
void SetParameters(const GlopParameters ¶meters)
Sets or gets the algorithm parameters to be used on the next Solve().
ABSL_MUST_USE_RESULT Status Solve(const LinearProgram &lp, TimeLimit *time_limit)
double DeterministicTime() const
ABSL_MUST_USE_RESULT Status MinimizeFromTransposedMatrixWithSlack(const DenseRow &objective, Fractional objective_scaling_factor, Fractional objective_offset, TimeLimit *time_limit)
void SetRandom(absl::BitGenRef random)
int64_t GetNumberOfIterations() const
void SetStartingVariableValuesForNextSolve(const DenseRow &values)
static Status OK()
Improves readability but identical to 0-arg constructor.
@ ERROR_LU
The LU factorization of the current basis couldn't be computed.
StrictITISpan< RowIndex, const Fractional > ConstView
void resize(IntType size)
void SetParameters(const GlopParameters ¶meters)
Sets the algorithm parameters.
const DenseRow & GetDenseRow() const
Fractional Get(ColIndex col) const
Getters for the variable values.
const DenseBitRow & GetCanDecreaseBitRow() const
const DenseBitRow & GetCanIncreaseBitRow() const
const DenseBitRow & GetNotBasicBitRow() const
const DenseBitRow & GetIsRelevantBitRow() const
constexpr ColIndex kInvalidCol(-1)
std::string GetVariableTypeString(VariableType variable_type)
Returns the string representation of the VariableType enum.
StrictITIVector< RowIndex, ColIndex > RowToColMapping
Fractional Square(Fractional f)
StrictITIVector< RowIndex, bool > DenseBooleanColumn
Column of booleans.
Permutation< ColIndex > ColumnPermutation
std::string Stringify(const Fractional x, bool fraction)
Fractional PreciseScalarProduct(const DenseRowOrColumn &u, const DenseRowOrColumn2 &v)
Bitset64< ColIndex > DenseBitRow
Row of bits.
const DenseRow & Transpose(const DenseColumn &col)
Similar comment as the other Transpose() implementation above.
std::string GetVariableStatusString(VariableStatus status)
Returns the string representation of the VariableStatus enum.
Index ColToIntIndex(ColIndex col)
Get the integer index corresponding to the col.
constexpr const uint64_t kDeterministicSeed
ColIndex RowToColIndex(RowIndex row)
Get the ColIndex corresponding to the column # row.
bool IsFinite(Fractional value)
constexpr RowIndex kInvalidRow(-1)
VariableType
Different types of variables.
@ UPPER_AND_LOWER_BOUNDED
bool AreFirstColumnsAndRowsExactlyEquals(RowIndex num_rows, ColIndex num_cols, const SparseMatrix &matrix_a, const CompactSparseMatrix &matrix_b)
StrictITIVector< ColIndex, VariableType > VariableTypeRow
Row of variable types.
Fractional InfinityNorm(const DenseColumn &v)
Returns the maximum of the coefficients of 'v'.
void ChangeSign(StrictITIVector< IndexType, Fractional > *data)
Changes the sign of all the entries in the given vector.
StrictITIVector< ColIndex, VariableStatus > VariableStatusRow
Row of variable statuses.
std::string StringifyMonomial(const Fractional a, absl::string_view x, bool fraction)
constexpr Fractional kInfinity
Infinity for type Fractional.
StrictITIVector< RowIndex, Fractional > DenseColumn
Column-vector types. Column-vector types are indexed by a row index.
StrictITIVector< ColIndex, Fractional > DenseRow
Row-vector types. Row-vector types are indexed by a column index.
ProblemStatus
Different statuses for a given problem.
@ ABNORMAL
An error occurred during the solving process.
@ INIT
The solver didn't had a chance to prove anything.
ConstraintStatus VariableToConstraintStatus(VariableStatus status)
Returns the ConstraintStatus corresponding to a given VariableStatus.
static double DeterministicTimeForFpOperations(int64_t n)
void ApplyColumnPermutationToRowIndexedVector(StrictITISpan< ColIndex, const ColIndex > col_perm, RowIndexedVector *v, RowIndexedVector *tmp)
std::string GetProblemStatusString(ProblemStatus problem_status)
Returns the string representation of the ProblemStatus enum.
In SWIG mode, we don't want anything besides these top-level includes.
util_intops::StrongVector< RowEntryIndex, SubsetIndex > SparseRow
DisabledScopedTimeDistributionUpdater ScopedTimeDistributionUpdater
#define RETURN_IF_NULL(x)
#define DCHECK_ROW_BOUNDS(row)
ABSL_FLAG(bool, simplex_display_numbers_as_fractions, false, "Display numbers as fractions.")
#define DCHECK_COL_BOUNDS(col)
#define IF_STATS_ENABLED(instructions)
#define SCOPED_TIME_STAT(stats)
#define GLOP_RETURN_IF_ERROR(function_call)
Macro to simplify error propagation between function returning Status.
#define GLOP_RETURN_ERROR_IF_NULL(arg)
Macro to check that a pointer argument is not null.
VariableStatusRow statuses
#define SOLVER_LOG(logger,...)