26#include "absl/flags/flag.h"
27#include "absl/log/check.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"
39#include "ortools/glop/parameters.pb.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); \
99bool UseAbslRandom() {
return false; }
110 random_(UseAbslRandom() ?
absl::BitGenRef(absl_random_)
111 :
absl::BitGenRef(deterministic_random_)),
112 basis_factorization_(&compact_matrix_, &basis_),
113 variables_info_(compact_matrix_),
114 primal_edge_norms_(compact_matrix_, variables_info_,
115 basis_factorization_),
116 dual_edge_norms_(basis_factorization_),
117 dual_prices_(random_),
118 variable_values_(parameters_, compact_matrix_, basis_, variables_info_,
119 basis_factorization_, &dual_edge_norms_, &dual_prices_),
120 update_row_(compact_matrix_, transposed_matrix_, variables_info_, basis_,
121 basis_factorization_),
122 reduced_costs_(compact_matrix_, objective_, basis_, variables_info_,
123 basis_factorization_, random_),
124 entering_variable_(variables_info_, random_, &reduced_costs_),
125 primal_prices_(random_, variables_info_, &primal_edge_norms_,
129 function_stats_(
"SimplexFunctionStats"),
137 solution_state_.statuses.clear();
138 variable_starting_values_.clear();
146 if (state.
statuses == solution_state_.statuses)
return;
148 solution_state_ = state;
149 solution_state_has_been_set_externally_ =
true;
154 variable_starting_values_ = values;
160 const double start_time =
time_limit->GetElapsedTime();
161 default_logger_.EnableLogging(parameters_.log_search_progress());
162 default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
163 parameters_ = initial_parameters_;
164 PropagateParameters();
167 if (transpose_was_changed_) {
168 compact_matrix_.PopulateFromTranspose(transposed_matrix_);
169 num_rows_ = compact_matrix_.num_rows();
170 num_cols_ = compact_matrix_.num_cols();
174 DCHECK_EQ(num_cols_, objective.
size());
177 objective_scaling_factor_ = objective_scaling_factor;
178 objective_offset_ = objective_offset;
179 const bool objective_is_unchanged = objective_ == objective;
180 objective_ = objective;
181 InitializeObjectiveLimit();
184 variables_info_.InitializeFromMutatedState();
186 if (objective_is_unchanged && parameters_.use_dual_simplex() &&
187 !transpose_was_changed_ && !solution_state_has_been_set_externally_ &&
188 !solution_state_.IsEmpty()) {
190 primal_edge_norms_.Clear();
191 variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
193 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
194 variable_values_.RecomputeBasicVariableValues();
195 return SolveInternal(start_time,
false, objective,
time_limit);
200 return SolveInternal(start_time,
false, objective,
time_limit);
204 const double start_time =
time_limit->GetElapsedTime();
205 default_logger_.EnableLogging(parameters_.log_search_progress());
206 default_logger_.SetLogToStdOut(parameters_.log_to_stdout());
214ABSL_MUST_USE_RESULT
Status RevisedSimplex::SolveInternal(
215 double start_time,
bool is_maximization_problem,
219 Cleanup update_deterministic_time_on_return(
227 DisplayBasicVariableStatistics();
230 dual_infeasibility_improvement_direction_.clear();
234 phase_ = Phase::FEASIBILITY;
236 num_feasibility_iterations_ = 0;
237 num_optimization_iterations_ = 0;
238 num_push_iterations_ = 0;
239 feasibility_time_ = 0.0;
240 optimization_time_ = 0.0;
248 solution_state_has_been_set_externally_ =
true;
251 ComputeNumberOfEmptyRows();
252 ComputeNumberOfEmptyColumns();
255 if (absl::GetFlag(FLAGS_simplex_stop_after_first_basis)) {
260 const bool use_dual = parameters_.use_dual_simplex();
265 primal_edge_norms_.SetPricingRule(parameters_.feasibility_rule());
267 if (parameters_.perturb_costs_in_dual_simplex()) {
268 reduced_costs_.PerturbCosts();
271 if (parameters_.use_dedicated_dual_feasibility_algorithm()) {
272 variables_info_.MakeBoxedVariableRelevant(
false);
274 DualMinimize(phase_ == Phase::FEASIBILITY,
time_limit));
287 variables_info_.MakeBoxedVariableRelevant(
true);
288 reduced_costs_.MakeReducedCostsPrecise();
291 MakeBoxedVariableDualFeasible(
292 variables_info_.GetNonBasicBoxedVariables(),
294 variable_values_.RecomputeBasicVariableValues();
300 reduced_costs_.MakeReducedCostsPrecise();
301 bool refactorize = reduced_costs_.NeedsBasisRefactorization();
305 reduced_costs_.ComputeMaximumDualInfeasibilityOnNonBoxedVariables();
306 if (initial_infeasibility <
307 reduced_costs_.GetDualFeasibilityTolerance()) {
308 SOLVER_LOG(logger_,
"Initial basis is dual feasible.");
310 MakeBoxedVariableDualFeasible(
311 variables_info_.GetNonBasicBoxedVariables(),
313 variable_values_.RecomputeBasicVariableValues();
316 variables_info_.TransformToDualPhaseIProblem(
317 reduced_costs_.GetDualFeasibilityTolerance(),
318 reduced_costs_.GetReducedCosts());
320 variable_values_.ResetAllNonBasicVariableValues(zero);
321 variable_values_.RecomputeBasicVariableValues();
329 variables_info_.EndDualPhaseI(
330 reduced_costs_.GetDualFeasibilityTolerance(),
331 reduced_costs_.GetFullReducedCosts());
332 variable_values_.ResetAllNonBasicVariableValues(
333 variable_starting_values_);
334 variable_values_.RecomputeBasicVariableValues();
344 if (reduced_costs_.ComputeMaximumDualInfeasibility() <
345 reduced_costs_.GetDualFeasibilityTolerance() + 1e-6) {
348 SOLVER_LOG(logger_,
"Infeasible after first phase.");
359 objective_ = objective_coefficients;
360 if (is_maximization_problem) {
365 objective_.
resize(num_cols_, 0.0);
366 reduced_costs_.ResetForNewObjective();
372 phase_ = Phase::OPTIMIZATION;
373 feasibility_time_ =
time_limit->GetElapsedTime() - start_time;
374 primal_edge_norms_.SetPricingRule(parameters_.optimization_rule());
375 num_feasibility_iterations_ = num_iterations_;
390 for (
int num_optims = 0;
394 num_optims <= parameters_.max_number_of_reoptimizations() &&
395 !objective_limit_reached_ &&
396 (num_iterations_ == 0 ||
397 num_iterations_ < parameters_.max_number_of_iterations()) &&
399 !absl::GetFlag(FLAGS_simplex_stop_after_feasibility) &&
409 DualMinimize(phase_ == Phase::FEASIBILITY,
time_limit));
417 basis_factorization_.IsRefactorized());
420 if (!integrality_scale_.empty() &&
429 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
432 variable_values_.RecomputeBasicVariableValues();
433 reduced_costs_.ClearAndRemoveCostShifts();
447 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
448 if (reduced_costs_.ComputeMaximumDualResidual() > tolerance ||
449 variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
450 variable_values_.ComputeMaximumPrimalInfeasibility() > tolerance) {
452 "PRIMAL_UNBOUNDED was reported, but the residual and/or "
453 "dual infeasibility is above the tolerance");
454 if (parameters_.change_status_to_imprecise()) {
475 double max_magnitude = 0.0;
477 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
478 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
479 double cost_delta = 0.0;
480 for (ColIndex col(0); col < num_cols_; ++col) {
481 cost_delta += solution_primal_ray_[col] * objective_[col];
482 if (solution_primal_ray_[col] > 0 && upper_bounds[col] !=
kInfinity) {
483 const Fractional value = variable_values_.Get(col);
484 const Fractional distance = (upper_bounds[col] - value + tolerance) /
485 solution_primal_ray_[col];
486 min_distance = std::min(distance, min_distance);
487 max_magnitude = std::max(solution_primal_ray_[col], max_magnitude);
489 if (solution_primal_ray_[col] < 0 && lower_bounds[col] != -
kInfinity) {
490 const Fractional value = variable_values_.Get(col);
491 const Fractional distance = (value - lower_bounds[col] + tolerance) /
492 -solution_primal_ray_[col];
493 min_distance = std::min(distance, min_distance);
494 max_magnitude = std::max(-solution_primal_ray_[col], max_magnitude);
497 SOLVER_LOG(logger_,
"Primal unbounded ray: max blocking magnitude = ",
498 max_magnitude,
", min distance to bound + ", tolerance,
" = ",
499 min_distance,
", ray cost delta = ", cost_delta);
500 if (min_distance * std::abs(cost_delta) < 1 &&
501 reduced_costs_.ComputeMaximumDualInfeasibility() <= tolerance) {
503 "PRIMAL_UNBOUNDED was reported, but the tolerance are good "
504 "and the unbounded ray is not great.");
506 "The difference between unbounded and optimal can depend "
507 "on a slight change of tolerance, trying to see if we are "
508 "at OPTIMAL after postsolve.");
514 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
515 if (reduced_costs_.ComputeMaximumDualResidual() > tolerance ||
516 variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
517 reduced_costs_.ComputeMaximumDualInfeasibility() > tolerance) {
519 "DUAL_UNBOUNDED was reported, but the residual and/or "
520 "dual infeasibility is above the tolerance");
521 if (parameters_.change_status_to_imprecise()) {
530 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
531 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
534 for (ColIndex col(0); col < num_cols_; ++col) {
535 const Fractional coeff = solution_dual_ray_row_combination_[col];
538 error = std::max(error, coeff);
540 implied_lb += coeff * lower_bounds[col];
542 }
else if (coeff < 0) {
544 error = std::max(error, -coeff);
546 implied_lb += coeff * upper_bounds[col];
551 " infeasibility=", implied_lb);
552 if (implied_lb < tolerance || error > tolerance) {
554 "DUAL_UNBOUNDED was reported, but the dual ray is not "
555 "proving infeasibility with high enough tolerance");
556 if (parameters_.change_status_to_imprecise()) {
567 parameters_.solution_feasibility_tolerance();
569 variable_values_.ComputeMaximumPrimalResidual();
571 reduced_costs_.ComputeMaximumDualResidual();
572 if (primal_residual > solution_tolerance ||
573 dual_residual > solution_tolerance) {
575 "OPTIMAL was reported, yet one of the residuals is "
576 "above the solution feasibility tolerance after the "
577 "shift/perturbation are removed.");
578 if (parameters_.change_status_to_imprecise()) {
588 primal_residual, parameters_.primal_feasibility_tolerance());
590 std::max(dual_residual, parameters_.dual_feasibility_tolerance());
592 variable_values_.ComputeMaximumPrimalInfeasibility();
594 reduced_costs_.ComputeMaximumDualInfeasibility();
595 if (primal_infeasibility > primal_tolerance &&
596 dual_infeasibility > dual_tolerance) {
598 "OPTIMAL was reported, yet both of the infeasibility "
599 "are above the tolerance after the "
600 "shift/perturbation are removed.");
601 if (parameters_.change_status_to_imprecise()) {
604 }
else if (primal_infeasibility > primal_tolerance) {
605 if (num_optims == parameters_.max_number_of_reoptimizations()) {
607 "The primal infeasibility is still higher than the "
608 "requested internal tolerance, but the maximum "
609 "number of optimization is reached.");
613 SOLVER_LOG(logger_,
"Re-optimizing with dual simplex ... ");
615 }
else if (dual_infeasibility > dual_tolerance) {
616 if (num_optims == parameters_.max_number_of_reoptimizations()) {
618 "The dual infeasibility is still higher than the "
619 "requested internal tolerance, but the maximum "
620 "number of optimization is reached.");
624 SOLVER_LOG(logger_,
"Re-optimizing with primal simplex ... ");
635 if (parameters_.change_status_to_imprecise() &&
637 const Fractional tolerance = parameters_.solution_feasibility_tolerance();
638 if (variable_values_.ComputeMaximumPrimalResidual() > tolerance ||
639 reduced_costs_.ComputeMaximumDualResidual() > tolerance) {
644 if (reduced_costs_.ComputeMaximumDualInfeasibility() > tolerance) {
650 if (variable_values_.ComputeMaximumPrimalInfeasibility() > tolerance) {
656 total_time_ =
time_limit->GetElapsedTime() - start_time;
657 optimization_time_ = total_time_ - feasibility_time_;
658 num_optimization_iterations_ = num_iterations_ - num_feasibility_iterations_;
662 if (!variable_starting_values_.empty()) {
663 const int num_super_basic = ComputeNumberOfSuperBasicVariables();
664 if (num_super_basic > 0) {
666 "Num super-basic variables left after optimize phase: ",
668 if (parameters_.push_to_vertex()) {
671 phase_ = Phase::PUSH;
677 "Skipping push phase because optimize didn't succeed.");
683 total_time_ =
time_limit->GetElapsedTime() - start_time;
684 push_time_ = total_time_ - feasibility_time_ - optimization_time_;
685 num_push_iterations_ = num_iterations_ - num_feasibility_iterations_ -
686 num_optimization_iterations_;
689 solution_objective_value_ = ComputeInitialProblemObjectiveValue();
690 solution_dual_values_ = reduced_costs_.GetDualValues();
691 solution_reduced_costs_ = reduced_costs_.GetReducedCosts();
694 if (is_maximization_problem) {
702 solution_objective_value_ =
705 if (is_maximization_problem) {
706 solution_objective_value_ = -solution_objective_value_;
710 variable_starting_values_.clear();
716 return problem_status_;
720 return solution_objective_value_;
724 return num_iterations_;
732 return variable_values_.Get(col);
736 return solution_reduced_costs_[col];
740 return solution_reduced_costs_;
744 return solution_dual_values_[row];
748 return variables_info_.GetStatusRow()[col];
756 return -variable_values_.Get(SlackColIndex(row));
762 const VariableStatus s = variables_info_.GetStatusRow()[SlackColIndex(row)];
774 return solution_primal_ray_;
778 return solution_dual_ray_;
783 return solution_dual_ray_row_combination_;
789 DCHECK(basis_factorization_.GetColumnPermutation().empty());
790 return basis_factorization_;
793std::string RevisedSimplex::GetPrettySolverStats()
const {
794 return absl::StrFormat(
795 "Problem status : %s\n"
796 "Solving time : %-6.4g\n"
797 "Number of iterations : %u\n"
798 "Time for solvability (first phase) : %-6.4g\n"
799 "Number of iterations for solvability : %u\n"
800 "Time for optimization : %-6.4g\n"
801 "Number of iterations for optimization : %u\n"
802 "Stop after first basis : %d\n",
804 feasibility_time_, num_feasibility_iterations_, optimization_time_,
805 num_optimization_iterations_,
806 absl::GetFlag(FLAGS_simplex_stop_after_first_basis));
812 basis_factorization_.DeterministicTime() +
813 update_row_.DeterministicTime() +
814 entering_variable_.DeterministicTime() +
815 reduced_costs_.DeterministicTime() +
816 primal_edge_norms_.DeterministicTime();
819void RevisedSimplex::SetVariableNames() {
820 variable_name_.
resize(num_cols_,
"");
821 for (ColIndex col(0); col < first_slack_col_; ++col) {
822 const ColIndex var_index = col + 1;
823 variable_name_[col] = absl::StrFormat(
"x%d",
ColToIntIndex(var_index));
825 for (ColIndex col(first_slack_col_); col < num_cols_; ++col) {
826 const ColIndex var_index = col - first_slack_col_ + 1;
827 variable_name_[col] = absl::StrFormat(
"s%d",
ColToIntIndex(var_index));
831void RevisedSimplex::SetNonBasicVariableStatusAndDeriveValue(
833 variables_info_.UpdateToNonBasicStatus(col, status);
834 variable_values_.SetNonBasicVariableValueFromStatus(col);
837bool RevisedSimplex::BasisIsConsistent()
const {
838 const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow();
840 for (RowIndex row(0); row < num_rows_; ++row) {
841 const ColIndex col = basis_[row];
842 if (!is_basic.IsSet(col))
return false;
845 ColIndex cols_in_basis(0);
846 ColIndex cols_not_in_basis(0);
847 for (ColIndex col(0); col < num_cols_; ++col) {
848 cols_in_basis += is_basic.IsSet(col);
849 cols_not_in_basis += !is_basic.IsSet(col);
850 if (is_basic.IsSet(col) !=
856 if (cols_not_in_basis != num_cols_ -
RowToColIndex(num_rows_))
return false;
862void RevisedSimplex::UpdateBasis(ColIndex entering_col, RowIndex basis_row,
870 DCHECK(!variables_info_.GetIsBasicBitRow().IsSet(entering_col));
871 DCHECK_NE(basis_[basis_row], entering_col);
874 const ColIndex leaving_col = basis_[basis_row];
875 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(leaving_col));
880 variables_info_.UpdateToNonBasicStatus(leaving_col, leaving_variable_status);
885 basis_[basis_row] = entering_col;
886 variables_info_.UpdateToBasicStatus(entering_col);
887 update_row_.Invalidate();
893class ColumnComparator {
895 explicit ColumnComparator(
const DenseRow& value) : value_(value) {}
896 bool operator()(ColIndex col_a, ColIndex col_b)
const {
897 return value_[col_a] < value_[col_b];
914void RevisedSimplex::UseSingletonColumnInInitialBasis(
RowToColMapping* basis) {
921 std::vector<ColIndex> singleton_column;
922 DenseRow cost_variation(num_cols_, 0.0);
923 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
924 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
925 for (ColIndex col(0); col < num_cols_; ++col) {
926 if (compact_matrix_.column(col).num_entries() != 1)
continue;
927 if (lower_bounds[col] == upper_bounds[col])
continue;
928 const Fractional slope = compact_matrix_.column(col).GetFirstCoefficient();
929 if (variable_values_.Get(col) == lower_bounds[col]) {
930 cost_variation[col] = objective_[col] / std::abs(slope);
932 cost_variation[col] = -objective_[col] / std::abs(slope);
934 singleton_column.push_back(col);
936 if (singleton_column.empty())
return;
943 ColumnComparator comparator(cost_variation);
944 std::sort(singleton_column.begin(), singleton_column.end(), comparator);
945 DCHECK_LE(cost_variation[singleton_column.front()],
946 cost_variation[singleton_column.back()]);
953 const DenseRow& variable_values = variable_values_.GetDenseRow();
954 for (
const ColIndex col : singleton_column) {
955 const RowIndex row = compact_matrix_.column(col).EntryRow(EntryIndex(0));
966 if (error_[row] == 0.0)
continue;
972 compact_matrix_.column(col).EntryCoefficient(EntryIndex(0));
973 const Fractional new_value = variable_values[col] + error_[row] / coeff;
974 if (new_value >= lower_bounds[col] && new_value <= upper_bounds[col]) {
985 const Fractional box_width = variables_info_.GetBoundDifference(col);
986 DCHECK_NE(box_width, 0.0);
987 DCHECK_NE(error_[row], 0.0);
988 const Fractional error_sign = error_[row] / coeff;
989 if (variable_values[col] == lower_bounds[col] && error_sign > 0.0) {
991 error_[row] -= coeff * box_width;
992 SetNonBasicVariableStatusAndDeriveValue(col,
996 if (variable_values[col] == upper_bounds[col] && error_sign < 0.0) {
998 error_[row] += coeff * box_width;
999 SetNonBasicVariableStatusAndDeriveValue(col,
1006bool RevisedSimplex::InitializeMatrixAndTestIfUnchanged(
1008 bool* only_change_is_new_rows,
bool* only_change_is_new_cols,
1009 ColIndex* num_new_cols) {
1011 DCHECK(only_change_is_new_rows !=
nullptr);
1012 DCHECK(only_change_is_new_cols !=
nullptr);
1013 DCHECK(num_new_cols !=
nullptr);
1014 DCHECK_EQ(num_cols_, compact_matrix_.num_cols());
1015 DCHECK_EQ(num_rows_, compact_matrix_.num_rows());
1018 const bool old_part_of_matrix_is_unchanged =
1020 num_rows_, first_slack_col_, lp.GetSparseMatrix(), compact_matrix_);
1023 const ColIndex lp_first_slack =
1024 lp_is_in_equation_form ? lp.GetFirstSlackVariable() : lp.num_variables();
1029 if (old_part_of_matrix_is_unchanged && lp.num_constraints() == num_rows_ &&
1030 lp_first_slack == first_slack_col_) {
1035 if (parameters_.use_transposed_matrix()) {
1036 if (transposed_matrix_.IsEmpty()) {
1037 transposed_matrix_.PopulateFromTranspose(compact_matrix_);
1040 transposed_matrix_.Reset(RowIndex(0));
1047 *only_change_is_new_rows = old_part_of_matrix_is_unchanged &&
1048 lp.num_constraints() > num_rows_ &&
1049 lp_first_slack == first_slack_col_;
1053 *only_change_is_new_cols = old_part_of_matrix_is_unchanged &&
1054 lp.num_constraints() == num_rows_ &&
1055 lp_first_slack > first_slack_col_;
1056 *num_new_cols = *only_change_is_new_cols ? lp_first_slack - first_slack_col_
1060 first_slack_col_ = lp_first_slack;
1063 num_rows_ = lp.num_constraints();
1064 num_cols_ = lp_first_slack +
RowToColIndex(lp.num_constraints());
1067 if (lp_is_in_equation_form) {
1070 compact_matrix_.PopulateFromMatrixView(MatrixView(lp.GetSparseMatrix()));
1072 compact_matrix_.PopulateFromSparseMatrixAndAddSlacks(lp.GetSparseMatrix());
1074 if (parameters_.use_transposed_matrix()) {
1075 transposed_matrix_.PopulateFromTranspose(compact_matrix_);
1077 transposed_matrix_.Reset(RowIndex(0));
1084bool RevisedSimplex::OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
1086 ColIndex num_new_cols) {
1088 DCHECK_LE(num_new_cols, first_slack_col_);
1089 const ColIndex first_new_col(first_slack_col_ - num_new_cols);
1092 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1093 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1094 for (ColIndex col(0); col < first_new_col; ++col) {
1095 if (lower_bounds[col] != lp.variable_lower_bounds()[col] ||
1096 upper_bounds[col] != lp.variable_upper_bounds()[col]) {
1102 for (ColIndex col(first_new_col); col < first_slack_col_; ++col) {
1103 if (lp.variable_lower_bounds()[col] != 0.0 &&
1104 lp.variable_upper_bounds()[col] != 0.0) {
1110 if (lp_is_in_equation_form) {
1111 for (ColIndex col(first_slack_col_); col < num_cols_; ++col) {
1112 if (lower_bounds[col - num_new_cols] != lp.variable_lower_bounds()[col] ||
1113 upper_bounds[col - num_new_cols] != lp.variable_upper_bounds()[col]) {
1118 DCHECK_EQ(num_rows_, lp.num_constraints());
1119 for (RowIndex row(0); row < num_rows_; ++row) {
1121 if (lower_bounds[col - num_new_cols] !=
1122 -lp.constraint_upper_bounds()[row] ||
1123 upper_bounds[col - num_new_cols] !=
1124 -lp.constraint_lower_bounds()[row]) {
1132bool RevisedSimplex::InitializeObjectiveAndTestIfUnchanged(
1136 bool objective_is_unchanged =
true;
1137 objective_.resize(num_cols_, 0.0);
1141 DCHECK_GE(num_cols_, lp.num_variables());
1143 const auto obj_coeffs = lp.objective_coefficients().const_view();
1144 for (ColIndex col(obj_coeffs.size()); col < num_cols_; ++col) {
1145 if (objective_[col] != 0.0) {
1146 objective_is_unchanged =
false;
1147 objective_[col] = 0.0;
1151 if (lp.IsMaximizationProblem()) {
1153 for (ColIndex col(0); col < obj_coeffs.size(); ++col) {
1155 if (objective_[col] != coeff) {
1156 objective_is_unchanged =
false;
1157 objective_[col] = coeff;
1160 objective_offset_ = -lp.objective_offset();
1161 objective_scaling_factor_ = -lp.objective_scaling_factor();
1163 for (ColIndex col(0); col < obj_coeffs.size(); ++col) {
1165 if (objective_[col] != coeff) {
1166 objective_is_unchanged =
false;
1167 objective_[col] = coeff;
1170 objective_offset_ = lp.objective_offset();
1171 objective_scaling_factor_ = lp.objective_scaling_factor();
1174 return objective_is_unchanged;
1177void RevisedSimplex::InitializeObjectiveLimit() {
1178 objective_limit_reached_ =
false;
1179 DCHECK(std::isfinite(objective_offset_));
1180 DCHECK(std::isfinite(objective_scaling_factor_));
1181 DCHECK_NE(0.0, objective_scaling_factor_);
1184 for (
const bool set_dual : {
true,
false}) {
1196 const Fractional limit = (objective_scaling_factor_ >= 0.0) != set_dual
1197 ? parameters_.objective_lower_limit()
1198 : parameters_.objective_upper_limit();
1200 limit / objective_scaling_factor_ - objective_offset_;
1202 dual_objective_limit_ = shifted_limit;
1204 primal_objective_limit_ = shifted_limit;
1214Status RevisedSimplex::CreateInitialBasis() {
1221 variables_info_.InitializeToDefaultStatus();
1222 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
1226 for (RowIndex row(0); row < num_rows_; ++row) {
1227 basis[row] = SlackColIndex(row);
1232 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1233 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1234 if (!parameters_.use_dual_simplex() &&
1235 parameters_.initial_basis() != GlopParameters::MAROS &&
1236 parameters_.exploit_singleton_column_in_initial_basis()) {
1240 for (ColIndex col(0); col < num_cols_; ++col) {
1241 if (compact_matrix_.column(col).num_entries() != 1)
continue;
1242 const VariableStatus status = variables_info_.GetStatusRow()[col];
1243 const Fractional objective = objective_[col];
1244 if (objective > 0 &&
IsFinite(lower_bounds[col]) &&
1246 SetNonBasicVariableStatusAndDeriveValue(col,
1248 }
else if (objective < 0 &&
IsFinite(upper_bounds[col]) &&
1250 SetNonBasicVariableStatusAndDeriveValue(col,
1257 ComputeVariableValuesError();
1266 UseSingletonColumnInInitialBasis(&basis);
1269 for (RowIndex row(0); row < num_rows_; ++row) {
1271 basis[row] = SlackColIndex(row);
1277 if (parameters_.initial_basis() == GlopParameters::NONE) {
1278 return InitializeFirstBasis(basis);
1280 if (parameters_.initial_basis() == GlopParameters::MAROS) {
1281 InitialBasis initial_basis(compact_matrix_, objective_, lower_bounds,
1282 upper_bounds, variables_info_.GetTypeRow());
1283 if (parameters_.use_dual_simplex()) {
1286 initial_basis.GetDualMarosBasis(num_cols_, &basis);
1288 initial_basis.GetPrimalMarosBasis(num_cols_, &basis);
1290 int number_changed = 0;
1291 for (RowIndex row(0); row < num_rows_; ++row) {
1292 if (basis[row] != SlackColIndex(row)) {
1296 VLOG(1) <<
"Number of Maros basis changes: " << number_changed;
1297 }
else if (parameters_.initial_basis() == GlopParameters::BIXBY ||
1298 parameters_.initial_basis() == GlopParameters::TRIANGULAR) {
1300 int num_fixed_variables = 0;
1301 for (RowIndex row(0); row < basis.size(); ++row) {
1302 const ColIndex col = basis[row];
1303 if (lower_bounds[col] == upper_bounds[col]) {
1305 ++num_fixed_variables;
1309 if (num_fixed_variables == 0) {
1310 SOLVER_LOG(logger_,
"Crash is set to ", parameters_.initial_basis(),
1311 " but there is no equality rows to remove from initial all "
1312 "slack basis. Starting from there.");
1315 SOLVER_LOG(logger_,
"Trying to remove ", num_fixed_variables,
1316 " fixed variables from the initial basis.");
1317 InitialBasis initial_basis(compact_matrix_, objective_, lower_bounds,
1318 upper_bounds, variables_info_.GetTypeRow());
1320 if (parameters_.initial_basis() == GlopParameters::BIXBY) {
1321 if (parameters_.use_scaling()) {
1322 initial_basis.CompleteBixbyBasis(first_slack_col_, &basis);
1324 VLOG(1) <<
"Bixby initial basis algorithm requires the problem "
1325 <<
"to be scaled. Skipping Bixby's algorithm.";
1327 }
else if (parameters_.initial_basis() == GlopParameters::TRIANGULAR) {
1330 if (parameters_.use_dual_simplex()) {
1333 initial_basis.CompleteTriangularDualBasis(num_cols_, &basis);
1335 initial_basis.CompleteTriangularPrimalBasis(num_cols_, &basis);
1338 const Status status = InitializeFirstBasis(basis);
1344 "Advanced basis algo failed, Reverting to all slack basis.");
1346 for (RowIndex row(0); row < num_rows_; ++row) {
1347 basis[row] = SlackColIndex(row);
1353 LOG(WARNING) <<
"Unsupported initial_basis parameters: "
1354 << parameters_.initial_basis();
1357 return InitializeFirstBasis(basis);
1366 for (RowIndex row(0); row < num_rows_; ++row) {
1368 basis_[row] = SlackColIndex(row);
1381 basis_factorization_.ComputeInfinityNormConditionNumberUpperBound();
1382 if (condition_number_ub > parameters_.initial_condition_number_threshold()) {
1383 const std::string error_message =
1384 absl::StrCat(
"The matrix condition number upper bound is too high: ",
1385 condition_number_ub);
1391 for (RowIndex row(0); row < num_rows_; ++row) {
1392 variables_info_.UpdateToBasicStatus(basis_[row]);
1394 DCHECK(BasisIsConsistent());
1396 variable_values_.ResetAllNonBasicVariableValues(variable_starting_values_);
1397 variable_values_.RecomputeBasicVariableValues();
1399 if (logger_->LoggingIsEnabled()) {
1403 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
1404 if (variable_values_.ComputeMaximumPrimalResidual() > tolerance) {
1407 "The primal residual of the initial basis is above the tolerance, ",
1408 variable_values_.ComputeMaximumPrimalResidual(),
" vs. ", tolerance);
1415 parameters_ = initial_parameters_;
1416 PropagateParameters();
1424 const bool lp_is_in_equation_form = lp.IsInEquationForm();
1431 ColIndex num_new_cols(0);
1432 bool only_change_is_new_rows =
false;
1433 bool only_change_is_new_cols =
false;
1434 const bool matrix_is_unchanged = InitializeMatrixAndTestIfUnchanged(
1435 lp, lp_is_in_equation_form, &only_change_is_new_rows,
1436 &only_change_is_new_cols, &num_new_cols);
1437 const bool only_new_bounds =
1438 only_change_is_new_cols && num_new_cols > 0 &&
1439 OldBoundsAreUnchangedAndNewVariablesHaveOneBoundAtZero(
1440 lp, lp_is_in_equation_form, num_new_cols);
1443 const bool objective_is_unchanged = InitializeObjectiveAndTestIfUnchanged(lp);
1445 const bool bounds_are_unchanged =
1446 lp_is_in_equation_form
1447 ? variables_info_.LoadBoundsAndReturnTrueIfUnchanged(
1448 lp.variable_lower_bounds(), lp.variable_upper_bounds())
1449 : variables_info_.LoadBoundsAndReturnTrueIfUnchanged(
1450 lp.variable_lower_bounds(), lp.variable_upper_bounds(),
1451 lp.constraint_lower_bounds(), lp.constraint_upper_bounds());
1456 if (matrix_is_unchanged && parameters_.allow_simplex_algorithm_change()) {
1457 if (objective_is_unchanged && !bounds_are_unchanged) {
1458 parameters_.set_use_dual_simplex(
true);
1459 PropagateParameters();
1461 if (bounds_are_unchanged && !objective_is_unchanged) {
1462 parameters_.set_use_dual_simplex(
false);
1463 PropagateParameters();
1467 InitializeObjectiveLimit();
1472 if (VLOG_IS_ON(2)) {
1483 bool solve_from_scratch =
true;
1486 if (!solution_state_.IsEmpty() && !solution_state_has_been_set_externally_) {
1487 if (!parameters_.use_dual_simplex()) {
1492 dual_edge_norms_.Clear();
1493 dual_pricing_vector_.clear();
1494 if (matrix_is_unchanged && bounds_are_unchanged) {
1497 reduced_costs_.ClearAndRemoveCostShifts();
1498 solve_from_scratch =
false;
1499 }
else if (only_change_is_new_cols && only_new_bounds) {
1500 variables_info_.InitializeFromBasisState(first_slack_col_, num_new_cols,
1502 variable_values_.ResetAllNonBasicVariableValues(
1503 variable_starting_values_);
1505 const ColIndex first_new_col(first_slack_col_ - num_new_cols);
1506 for (ColIndex& col_ref : basis_) {
1507 if (col_ref >= first_new_col) {
1508 col_ref += num_new_cols;
1515 primal_edge_norms_.Clear();
1516 reduced_costs_.ClearAndRemoveCostShifts();
1517 solve_from_scratch =
false;
1523 primal_edge_norms_.Clear();
1524 if (objective_is_unchanged) {
1525 if (matrix_is_unchanged) {
1526 if (!bounds_are_unchanged) {
1527 variables_info_.InitializeFromBasisState(
1528 first_slack_col_, ColIndex(0), solution_state_);
1529 variable_values_.ResetAllNonBasicVariableValues(
1530 variable_starting_values_);
1531 variable_values_.RecomputeBasicVariableValues();
1533 solve_from_scratch =
false;
1534 }
else if (only_change_is_new_rows) {
1537 variables_info_.InitializeFromBasisState(
1538 first_slack_col_, ColIndex(0), solution_state_);
1539 dual_edge_norms_.ResizeOnNewRows(num_rows_);
1544 reduced_costs_.ClearAndRemoveCostShifts();
1545 dual_pricing_vector_.clear();
1548 if (InitializeFirstBasis(basis_).ok()) {
1549 solve_from_scratch =
false;
1556 return FinishInitialization(solve_from_scratch);
1559Status RevisedSimplex::FinishInitialization(
bool solve_from_scratch) {
1562 if (solve_from_scratch && !solution_state_.IsEmpty()) {
1563 basis_factorization_.Clear();
1564 reduced_costs_.ClearAndRemoveCostShifts();
1565 primal_edge_norms_.Clear();
1566 dual_edge_norms_.Clear();
1567 dual_pricing_vector_.clear();
1571 variables_info_.InitializeFromBasisState(first_slack_col_, ColIndex(0),
1575 std::vector<ColIndex> candidates;
1576 for (
const ColIndex col : variables_info_.GetIsBasicBitRow()) {
1577 candidates.push_back(col);
1579 SOLVER_LOG(logger_,
"The warm-start state contains ", candidates.size(),
1580 " candidates for the basis (num_rows = ", num_rows_.value(),
1586 if (candidates.size() == num_rows_) {
1588 for (
const ColIndex col : candidates) {
1589 basis_.push_back(col);
1595 if (InitializeFirstBasis(basis_).ok()) {
1596 solve_from_scratch =
false;
1600 if (solve_from_scratch) {
1601 basis_ = basis_factorization_.ComputeInitialBasis(candidates);
1602 const int num_super_basic =
1603 variables_info_.ChangeUnusedBasicVariablesToFree(basis_);
1604 const int num_snapped = variables_info_.SnapFreeVariablesToBound(
1605 parameters_.crossover_bound_snapping_distance(),
1606 variable_starting_values_);
1607 if (logger_->LoggingIsEnabled()) {
1608 SOLVER_LOG(logger_,
"The initial basis did not use ",
1609 " BASIC columns from the initial state and used ",
1610 (num_rows_ - (candidates.size() - num_super_basic)).value(),
1611 " slack variables that were not marked BASIC.");
1612 if (num_snapped > 0) {
1614 " of the FREE variables where moved to their bound.");
1618 if (InitializeFirstBasis(basis_).ok()) {
1619 solve_from_scratch =
false;
1622 "RevisedSimplex is not using the warm start "
1623 "basis because it is not factorizable.");
1628 if (solve_from_scratch) {
1629 SOLVER_LOG(logger_,
"Starting basis: create from scratch.");
1630 basis_factorization_.Clear();
1631 reduced_costs_.ClearAndRemoveCostShifts();
1632 primal_edge_norms_.Clear();
1633 dual_edge_norms_.Clear();
1634 dual_pricing_vector_.clear();
1637 SOLVER_LOG(logger_,
"Starting basis: incremental solve.");
1639 DCHECK(BasisIsConsistent());
1641 transpose_was_changed_ =
false;
1645void RevisedSimplex::DisplayBasicVariableStatistics() {
1648 int num_fixed_variables = 0;
1649 int num_free_variables = 0;
1650 int num_variables_at_bound = 0;
1651 int num_slack_variables = 0;
1652 int num_infeasible_variables = 0;
1654 const DenseRow& variable_values = variable_values_.GetDenseRow();
1656 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1657 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1658 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
1659 for (RowIndex row(0); row < num_rows_; ++row) {
1660 const ColIndex col = basis_[row];
1661 const Fractional value = variable_values[col];
1663 ++num_free_variables;
1665 if (value > upper_bounds[col] + tolerance ||
1666 value < lower_bounds[col] - tolerance) {
1667 ++num_infeasible_variables;
1669 if (col >= first_slack_col_) {
1670 ++num_slack_variables;
1672 if (lower_bounds[col] == upper_bounds[col]) {
1673 ++num_fixed_variables;
1674 }
else if (variable_values[col] == lower_bounds[col] ||
1675 variable_values[col] == upper_bounds[col]) {
1676 ++num_variables_at_bound;
1680 SOLVER_LOG(logger_,
"The matrix with slacks has ",
1681 compact_matrix_.num_rows().value(),
" rows, ",
1682 compact_matrix_.num_cols().value(),
" columns, ",
1683 compact_matrix_.num_entries().value(),
" entries.");
1684 SOLVER_LOG(logger_,
"Number of basic infeasible variables: ",
1685 num_infeasible_variables);
1686 SOLVER_LOG(logger_,
"Number of basic slack variables: ", num_slack_variables);
1688 "Number of basic variables at bound: ", num_variables_at_bound);
1689 SOLVER_LOG(logger_,
"Number of basic fixed variables: ", num_fixed_variables);
1690 SOLVER_LOG(logger_,
"Number of basic free variables: ", num_free_variables);
1691 SOLVER_LOG(logger_,
"Number of super-basic variables: ",
1692 ComputeNumberOfSuperBasicVariables());
1695void RevisedSimplex::SaveState() {
1696 DCHECK_EQ(num_cols_, variables_info_.GetStatusRow().size());
1697 solution_state_.statuses = variables_info_.GetStatusRow();
1698 solution_state_has_been_set_externally_ =
false;
1701RowIndex RevisedSimplex::ComputeNumberOfEmptyRows() {
1703 for (ColIndex col(0); col < num_cols_; ++col) {
1704 for (
const SparseColumn::Entry e : compact_matrix_.column(col)) {
1705 contains_data[e.row()] =
true;
1708 RowIndex num_empty_rows(0);
1709 for (RowIndex row(0); row < num_rows_; ++row) {
1710 if (!contains_data[row]) {
1712 VLOG(2) <<
"Row " << row <<
" is empty.";
1715 return num_empty_rows;
1718ColIndex RevisedSimplex::ComputeNumberOfEmptyColumns() {
1719 ColIndex num_empty_cols(0);
1720 for (ColIndex col(0); col < num_cols_; ++col) {
1721 if (compact_matrix_.column(col).IsEmpty()) {
1723 VLOG(2) <<
"Column " << col <<
" is empty.";
1726 return num_empty_cols;
1729int RevisedSimplex::ComputeNumberOfSuperBasicVariables()
const {
1731 int num_super_basic = 0;
1732 for (ColIndex col(0); col < num_cols_; ++col) {
1734 variable_values_.Get(col) != 0.0) {
1738 return num_super_basic;
1741void RevisedSimplex::CorrectErrorsOnVariableValues() {
1743 DCHECK(basis_factorization_.IsRefactorized());
1749 variable_values_.ComputeMaximumPrimalResidual();
1753 if (primal_residual >= parameters_.harris_tolerance_ratio() *
1754 parameters_.primal_feasibility_tolerance()) {
1755 variable_values_.RecomputeBasicVariableValues();
1756 VLOG(1) <<
"Primal infeasibility (bounds error) = "
1757 << variable_values_.ComputeMaximumPrimalInfeasibility()
1758 <<
", Primal residual |A.x - b| = "
1759 << variable_values_.ComputeMaximumPrimalResidual();
1763void RevisedSimplex::ComputeVariableValuesError() {
1765 error_.AssignToZero(num_rows_);
1766 const DenseRow& variable_values = variable_values_.GetDenseRow();
1767 for (ColIndex col(0); col < num_cols_; ++col) {
1768 const Fractional value = variable_values[col];
1769 compact_matrix_.ColumnAddMultipleToDenseColumn(col, -value, &error_);
1773void RevisedSimplex::ComputeDirection(ColIndex col) {
1776 basis_factorization_.RightSolveForProblemColumn(col, &direction_);
1778 if (direction_.non_zeros.empty()) {
1780 const RowIndex num_rows = num_rows_;
1782 for (RowIndex row(0); row < num_rows; ++row) {
1785 direction_.non_zeros.push_back(row);
1786 norm = std::max(norm, std::abs(value));
1790 for (
const auto e : direction_) {
1791 norm = std::max(norm, std::abs(e.coefficient()));
1794 direction_infinity_norm_ = norm;
1796 num_rows_ == 0 ? 0.0
1797 :
static_cast<double>(direction_.non_zeros.size()) /
1798 static_cast<double>(num_rows_.value())));
1801Fractional RevisedSimplex::ComputeDirectionError(ColIndex col) {
1803 compact_matrix_.ColumnCopyToDenseColumn(col, &error_);
1804 for (
const auto e : direction_) {
1805 compact_matrix_.ColumnAddMultipleToDenseColumn(col, -e.coefficient(),
1811template <
bool is_entering_reduced_cost_positive>
1814 RowIndex row)
const {
1815 const ColIndex col = basis_[row];
1816 const Fractional direction = direction_[row];
1817 const Fractional value = variable_values_.Get(col);
1818 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(col));
1819 DCHECK_NE(direction, 0.0);
1820 if (is_entering_reduced_cost_positive) {
1821 if (direction > 0.0) {
1822 return (upper_bounds[col] - value) / direction;
1824 return (lower_bounds[col] - value) / direction;
1827 if (direction > 0.0) {
1828 return (value - lower_bounds[col]) / direction;
1830 return (value - upper_bounds[col]) / direction;
1835template <
bool is_entering_reduced_cost_positive>
1836Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates(
1840 parameters_.harris_tolerance_ratio() *
1841 parameters_.primal_feasibility_tolerance();
1842 const Fractional minimum_delta = parameters_.degenerate_ministep_factor() *
1843 parameters_.primal_feasibility_tolerance();
1849 leaving_candidates->Clear();
1855 const Fractional threshold = basis_factorization_.IsRefactorized()
1856 ? parameters_.minimum_acceptable_pivot()
1857 : parameters_.ratio_test_zero_threshold();
1859 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1860 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1861 for (
const auto e : direction_) {
1862 const Fractional magnitude = std::abs(e.coefficient());
1863 if (magnitude <= threshold)
continue;
1864 const Fractional ratio = GetRatio<is_entering_reduced_cost_positive>(
1865 lower_bounds, upper_bounds, e.row());
1866 if (ratio <= harris_ratio) {
1867 leaving_candidates->SetCoefficient(e.row(), ratio);
1879 harris_ratio = std::min(harris_ratio,
1880 std::max(minimum_delta / magnitude,
1881 ratio + harris_tolerance / magnitude));
1884 return harris_ratio;
1897 if (current >= 0.0) {
1898 return candidate >= 0.0 && candidate <= current;
1900 return candidate >= current;
1908Status RevisedSimplex::ChooseLeavingVariableRow(
1909 ColIndex entering_col,
Fractional reduced_cost,
bool* refactorize,
1916 DCHECK_NE(0.0, reduced_cost);
1919 int stats_num_leaving_choices = 0;
1920 equivalent_leaving_choices_.clear();
1921 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
1922 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
1924 stats_num_leaving_choices = 0;
1928 const Fractional entering_value = variable_values_.Get(entering_col);
1930 (reduced_cost > 0.0) ? entering_value - lower_bounds[entering_col]
1931 : upper_bounds[entering_col] - entering_value;
1932 DCHECK_GT(current_ratio, 0.0);
1938 (reduced_cost > 0.0) ? ComputeHarrisRatioAndLeavingCandidates<true>(
1939 current_ratio, &leaving_candidates_)
1940 : ComputeHarrisRatioAndLeavingCandidates<false>(
1941 current_ratio, &leaving_candidates_);
1946 if (current_ratio <= harris_ratio) {
1948 *step_length = current_ratio;
1958 stats_num_leaving_choices = 0;
1960 equivalent_leaving_choices_.clear();
1961 for (
const SparseColumn::Entry e : leaving_candidates_) {
1963 if (ratio > harris_ratio)
continue;
1964 ++stats_num_leaving_choices;
1965 const RowIndex row = e.row();
1970 const Fractional candidate_magnitude = std::abs(direction_[row]);
1971 if (candidate_magnitude < pivot_magnitude)
continue;
1972 if (candidate_magnitude == pivot_magnitude) {
1973 if (!IsRatioMoreOrEquallyStable(ratio, current_ratio))
continue;
1974 if (ratio == current_ratio) {
1976 equivalent_leaving_choices_.push_back(row);
1980 equivalent_leaving_choices_.clear();
1981 current_ratio = ratio;
1982 pivot_magnitude = candidate_magnitude;
1987 if (!equivalent_leaving_choices_.empty()) {
1988 equivalent_leaving_choices_.push_back(*leaving_row);
1990 equivalent_leaving_choices_[std::uniform_int_distribution<int>(
1991 0, equivalent_leaving_choices_.size() - 1)(random_)];
2003 if (current_ratio <= 0.0) {
2007 parameters_.degenerate_ministep_factor() *
2008 parameters_.primal_feasibility_tolerance();
2009 *step_length = minimum_delta / pivot_magnitude;
2011 *step_length = current_ratio;
2018 TestPivot(entering_col, *leaving_row);
2031 if (pivot_magnitude <
2032 parameters_.small_pivot_threshold() * direction_infinity_norm_) {
2036 if (!basis_factorization_.IsRefactorized()) {
2037 VLOG(1) <<
"Refactorizing to avoid pivoting by "
2038 << direction_[*leaving_row]
2039 <<
" direction_infinity_norm_ = " << direction_infinity_norm_
2040 <<
" reduced cost = " << reduced_cost;
2041 *refactorize =
true;
2051 VLOG(1) <<
"Couldn't avoid pivoting by " << direction_[*leaving_row]
2052 <<
" direction_infinity_norm_ = " << direction_infinity_norm_
2053 <<
" reduced cost = " << reduced_cost;
2054 DCHECK_GE(std::abs(direction_[*leaving_row]),
2055 parameters_.minimum_acceptable_pivot());
2063 const bool is_reduced_cost_positive = (reduced_cost > 0.0);
2064 const bool is_leaving_coeff_positive = (direction_[*leaving_row] > 0.0);
2065 *target_bound = (is_reduced_cost_positive == is_leaving_coeff_positive)
2066 ? upper_bounds[basis_[*leaving_row]]
2067 : lower_bounds[basis_[*leaving_row]];
2072 ratio_test_stats_.leaving_choices.Add(stats_num_leaving_choices);
2073 if (!equivalent_leaving_choices_.empty()) {
2074 ratio_test_stats_.num_perfect_ties.Add(
2075 equivalent_leaving_choices_.size());
2078 ratio_test_stats_.abs_used_pivot.Add(std::abs(direction_[*leaving_row]));
2094 coeff_magnitude(_coeff_magnitude),
2095 target_bound(_target_bound) {}
2100 bool operator<(
const BreakPoint& other)
const {
2101 if (ratio == other.ratio) {
2102 if (coeff_magnitude == other.coeff_magnitude) {
2103 return row > other.row;
2105 return coeff_magnitude < other.coeff_magnitude;
2107 return ratio > other.ratio;
2118void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow(
2119 ColIndex entering_col,
Fractional reduced_cost,
bool* refactorize,
2120 RowIndex* leaving_row,
Fractional* step_length,
2127 DCHECK_NE(0.0, reduced_cost);
2128 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2129 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2133 const Fractional entering_value = variable_values_.Get(entering_col);
2134 Fractional current_ratio = (reduced_cost > 0.0)
2135 ? entering_value - lower_bounds[entering_col]
2136 : upper_bounds[entering_col] - entering_value;
2137 DCHECK_GT(current_ratio, 0.0);
2139 std::vector<BreakPoint> breakpoints;
2140 const Fractional tolerance = parameters_.primal_feasibility_tolerance();
2141 for (
const auto e : direction_) {
2143 reduced_cost > 0.0 ? e.coefficient() : -e.coefficient();
2144 const Fractional magnitude = std::abs(direction);
2145 if (magnitude < tolerance)
continue;
2160 const ColIndex col = basis_[e.row()];
2161 DCHECK(variables_info_.GetIsBasicBitRow().IsSet(col));
2163 const Fractional value = variable_values_.Get(col);
2164 const Fractional lower_bound = lower_bounds[col];
2165 const Fractional upper_bound = upper_bounds[col];
2166 const Fractional to_lower = (lower_bound - tolerance - value) / direction;
2167 const Fractional to_upper = (upper_bound + tolerance - value) / direction;
2171 if (to_lower >= 0.0 && to_lower < current_ratio) {
2172 breakpoints.push_back(
2173 BreakPoint(e.row(), to_lower, magnitude, lower_bound));
2175 if (to_upper >= 0.0 && to_upper < current_ratio) {
2176 breakpoints.push_back(
2177 BreakPoint(e.row(), to_upper, magnitude, upper_bound));
2183 std::make_heap(breakpoints.begin(), breakpoints.end());
2187 Fractional improvement = std::abs(reduced_cost);
2190 while (!breakpoints.empty()) {
2191 const BreakPoint top = breakpoints.front();
2199 if (top.coeff_magnitude > best_magnitude) {
2200 *leaving_row = top.row;
2201 current_ratio = top.ratio;
2202 best_magnitude = top.coeff_magnitude;
2203 *target_bound = top.target_bound;
2208 improvement -= top.coeff_magnitude;
2209 if (improvement <= 0.0)
break;
2210 std::pop_heap(breakpoints.begin(), breakpoints.end());
2211 breakpoints.pop_back();
2217 parameters_.small_pivot_threshold() * direction_infinity_norm_;
2218 if (best_magnitude < threshold && !basis_factorization_.IsRefactorized()) {
2219 *refactorize =
true;
2223 *step_length = current_ratio;
2227Status RevisedSimplex::DualChooseLeavingVariableRow(RowIndex* leaving_row,
2236 if (dual_prices_.Size() == 0) {
2237 variable_values_.RecomputeDualPrices(
2238 parameters_.dual_price_prioritize_norm());
2243 *leaving_row = dual_prices_.GetMaximum();
2246 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2247 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2248 const ColIndex leaving_col = basis_[*leaving_row];
2249 const Fractional value = variable_values_.Get(leaving_col);
2250 if (value < lower_bounds[leaving_col]) {
2251 *cost_variation = lower_bounds[leaving_col] - value;
2252 *target_bound = lower_bounds[leaving_col];
2253 DCHECK_GT(*cost_variation, 0.0);
2255 *cost_variation = upper_bounds[leaving_col] - value;
2256 *target_bound = upper_bounds[leaving_col];
2257 DCHECK_LT(*cost_variation, 0.0);
2268 if (cost == 0.0)
return false;
2278template <
bool use_dense_update>
2282 const Fractional price = dual_pricing_vector_[row];
2283 const bool is_candidate =
2284 IsDualPhaseILeavingCandidate(price, type, threshold);
2286 if (use_dense_update) {
2287 dual_prices_.DenseAddOrUpdate(row,
Square(price) / squared_norm[row]);
2289 dual_prices_.AddOrUpdate(row,
Square(price) / squared_norm[row]);
2292 dual_prices_.Remove(row);
2296void RevisedSimplex::DualPhaseIUpdatePrice(RowIndex leaving_row,
2297 ColIndex entering_col) {
2306 if (reduced_costs_.AreReducedCostsRecomputed() ||
2307 dual_edge_norms_.NeedsBasisRefactorization() ||
2308 dual_pricing_vector_.empty()) {
2313 const Fractional threshold = parameters_.ratio_test_zero_threshold();
2318 dual_edge_norms_.GetEdgeSquaredNorms();
2324 dual_pricing_vector_[leaving_row] / direction_[leaving_row];
2325 for (
const auto e : direction_) {
2326 dual_pricing_vector_[e.row()] -= e.coefficient() * step;
2327 OnDualPriceChange(squared_norms, e.row(), variable_type[basis_[e.row()]],
2330 dual_pricing_vector_[leaving_row] = step;
2334 dual_pricing_vector_[leaving_row] -=
2335 dual_infeasibility_improvement_direction_[entering_col];
2336 if (dual_infeasibility_improvement_direction_[entering_col] != 0.0) {
2337 --num_dual_infeasible_positions_;
2339 dual_infeasibility_improvement_direction_[entering_col] = 0.0;
2342 dual_infeasibility_improvement_direction_[basis_[leaving_row]] = 0.0;
2345 OnDualPriceChange(squared_norms, leaving_row, variable_type[entering_col],
2349template <
typename Cols>
2350void RevisedSimplex::DualPhaseIUpdatePriceOnReducedCostChange(
2353 bool something_to_do =
false;
2354 const DenseBitRow::ConstView can_decrease =
2355 variables_info_.GetCanDecreaseBitRow().const_view();
2356 const DenseBitRow::ConstView can_increase =
2357 variables_info_.GetCanIncreaseBitRow().const_view();
2359 const Fractional tolerance = reduced_costs_.GetDualFeasibilityTolerance();
2360 auto improvement_direction = dual_infeasibility_improvement_direction_.view();
2361 for (
const ColIndex col : cols) {
2362 const Fractional reduced_cost = reduced_costs[col];
2364 (can_increase[col] && reduced_cost < -tolerance) ? 1.0
2365 : (can_decrease[col] && reduced_cost > tolerance) ? -1.0
2367 if (sign != improvement_direction[col]) {
2369 --num_dual_infeasible_positions_;
2370 }
else if (improvement_direction[col] == 0.0) {
2371 ++num_dual_infeasible_positions_;
2373 if (!something_to_do) {
2374 initially_all_zero_scratchpad_.values.resize(num_rows_, 0.0);
2375 initially_all_zero_scratchpad_.ClearSparseMask();
2376 initially_all_zero_scratchpad_.non_zeros.clear();
2377 something_to_do =
true;
2381 num_update_price_operations_ +=
2382 10 * compact_matrix_.column(col).num_entries().value();
2383 compact_matrix_.ColumnAddMultipleToSparseScatteredColumn(
2384 col, sign - improvement_direction[col],
2385 &initially_all_zero_scratchpad_);
2386 improvement_direction[col] = sign;
2389 if (something_to_do) {
2390 initially_all_zero_scratchpad_.ClearNonZerosIfTooDense();
2391 initially_all_zero_scratchpad_.ClearSparseMask();
2393 dual_edge_norms_.GetEdgeSquaredNorms();
2396 const Fractional threshold = parameters_.ratio_test_zero_threshold();
2397 basis_factorization_.RightSolve(&initially_all_zero_scratchpad_);
2398 if (initially_all_zero_scratchpad_.non_zeros.empty()) {
2399 dual_prices_.StartDenseUpdates();
2400 for (RowIndex row(0); row < num_rows_; ++row) {
2401 if (initially_all_zero_scratchpad_[row] == 0.0)
continue;
2402 dual_pricing_vector_[row] += initially_all_zero_scratchpad_[row];
2403 OnDualPriceChange<
true>(
2404 squared_norms, row, variable_type[basis_[row]], threshold);
2406 initially_all_zero_scratchpad_.values.AssignToZero(num_rows_);
2408 for (
const auto e : initially_all_zero_scratchpad_) {
2409 dual_pricing_vector_[e.row()] += e.coefficient();
2410 OnDualPriceChange(squared_norms, e.row(),
2411 variable_type[basis_[e.row()]], threshold);
2412 initially_all_zero_scratchpad_[e.row()] = 0.0;
2415 initially_all_zero_scratchpad_.non_zeros.clear();
2419Status RevisedSimplex::DualPhaseIChooseLeavingVariableRow(
2420 RowIndex* leaving_row,
Fractional* cost_variation,
2425 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2426 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2438 if (reduced_costs_.AreReducedCostsRecomputed() ||
2439 dual_edge_norms_.NeedsBasisRefactorization() ||
2440 dual_pricing_vector_.empty()) {
2442 num_dual_infeasible_positions_ = 0;
2443 dual_pricing_vector_.AssignToZero(num_rows_);
2444 dual_prices_.ClearAndResize(num_rows_);
2445 dual_infeasibility_improvement_direction_.AssignToZero(num_cols_);
2446 DualPhaseIUpdatePriceOnReducedCostChange(
2447 variables_info_.GetIsRelevantBitRow());
2451 DualPhaseIUpdatePriceOnReducedCostChange(update_row_.GetNonZeroPositions());
2456 if (num_dual_infeasible_positions_ == 0)
return Status::OK();
2458 *leaving_row = dual_prices_.GetMaximum();
2463 *cost_variation = dual_pricing_vector_[*leaving_row];
2464 const ColIndex leaving_col = basis_[*leaving_row];
2465 if (*cost_variation < 0.0) {
2466 *target_bound = upper_bounds[leaving_col];
2468 *target_bound = lower_bounds[leaving_col];
2474template <
typename BoxedVariableCols>
2475void RevisedSimplex::MakeBoxedVariableDualFeasible(
2476 const BoxedVariableCols& cols,
bool update_basic_values) {
2478 std::vector<ColIndex> changed_cols;
2488 const Fractional threshold = reduced_costs_.GetDualFeasibilityTolerance();
2490 const DenseRow& variable_values = variable_values_.GetDenseRow();
2492 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2493 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2495 for (
const ColIndex col : cols) {
2496 const Fractional reduced_cost = reduced_costs[col];
2498 DCHECK(variables_info_.GetTypeRow()[col] ==
2501 DCHECK(variable_values[col] == lower_bounds[col] ||
2502 variable_values[col] == upper_bounds[col] ||
2505 variables_info_.UpdateToNonBasicStatus(col,
2507 changed_cols.push_back(col);
2508 }
else if (reduced_cost < -threshold &&
2510 variables_info_.UpdateToNonBasicStatus(col,
2512 changed_cols.push_back(col);
2516 if (!changed_cols.empty()) {
2517 iteration_stats_.num_dual_flips.Add(changed_cols.size());
2518 variable_values_.UpdateGivenNonBasicVariables(changed_cols,
2519 update_basic_values);
2523Fractional RevisedSimplex::ComputeStepToMoveBasicVariableToBound(
2524 RowIndex leaving_row,
Fractional target_bound) {
2528 const ColIndex leaving_col = basis_[leaving_row];
2529 const Fractional leaving_variable_value = variable_values_.Get(leaving_col);
2530 Fractional unscaled_step = leaving_variable_value - target_bound;
2539 return unscaled_step / direction_[leaving_row];
2542bool RevisedSimplex::TestPivot(ColIndex entering_col, RowIndex leaving_row) {
2543 VLOG(1) <<
"Test pivot.";
2545 const ColIndex leaving_col = basis_[leaving_row];
2546 basis_[leaving_row] = entering_col;
2550 CompactSparseMatrixView basis_matrix(&compact_matrix_, &basis_);
2551 const bool is_ok = test_lu_.ComputeFactorization(basis_matrix).ok();
2552 basis_[leaving_row] = leaving_col;
2559void RevisedSimplex::PermuteBasis() {
2565 basis_factorization_.GetColumnPermutation();
2566 if (col_perm.empty())
return;
2573 if (!dual_pricing_vector_.empty()) {
2577 &dual_pricing_vector_,
2578 &tmp_dual_pricing_vector_);
2582 reduced_costs_.UpdateDataOnBasisPermutation();
2583 dual_edge_norms_.UpdateDataOnBasisPermutation(col_perm);
2587 basis_factorization_.SetColumnPermutationToIdentity();
2590Status RevisedSimplex::UpdateAndPivot(ColIndex entering_col,
2591 RowIndex leaving_row,
2607 if (update_row_.IsComputedFor(leaving_row)) {
2608 pivot_from_update_row = update_row_.GetCoefficient(entering_col);
2612 update_row_.ComputeUnitRowLeftInverse(leaving_row);
2613 pivot_from_update_row = compact_matrix_.ColumnScalarProduct(
2614 entering_col, update_row_.GetUnitRowLeftInverse().values);
2617 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
2618 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
2619 const ColIndex leaving_col = basis_[leaving_row];
2621 lower_bounds[leaving_col] == upper_bounds[leaving_col]
2623 : target_bound == lower_bounds[leaving_col]
2626 if (variable_values_.Get(leaving_col) != target_bound) {
2627 ratio_test_stats_.bound_shift.Add(variable_values_.Get(leaving_col) -
2630 UpdateBasis(entering_col, leaving_row, leaving_variable_status);
2633 const Fractional pivot_from_direction = direction_[leaving_row];
2635 std::abs(pivot_from_update_row - pivot_from_direction);
2636 if (diff > parameters_.refactorization_threshold() *
2637 (1.0 + std::min(std::abs(pivot_from_update_row),
2638 std::abs(pivot_from_direction)))) {
2639 VLOG(1) <<
"Refactorizing: imprecise pivot " << pivot_from_direction
2640 <<
" diff = " << diff;
2643 if (basis_factorization_.NumUpdates() < 10) {
2644 Fractional threshold = parameters_.lu_factorization_pivot_threshold();
2645 threshold = std::min(threshold * 1.5, 0.9);
2646 VLOG(1) <<
"Increasing LU pivot threshold " << threshold;
2647 parameters_.set_lu_factorization_pivot_threshold(threshold);
2648 basis_factorization_.SetParameters(parameters_);
2651 last_refactorization_reason_ = RefactorizationReason::IMPRECISE_PIVOT;
2655 basis_factorization_.Update(entering_col, leaving_row, direction_));
2657 if (basis_factorization_.IsRefactorized()) {
2663Status RevisedSimplex::RefactorizeBasisIfNeeded(
bool* refactorize) {
2665 if (*refactorize && !basis_factorization_.IsRefactorized()) {
2667 update_row_.Invalidate();
2670 *refactorize =
false;
2675 if (col >= integrality_scale_.size()) {
2676 integrality_scale_.resize(col + 1, 0.0);
2678 integrality_scale_[col] = scale;
2683 Cleanup update_deterministic_time_on_return(
2690 std::vector<ColIndex> candidates;
2693 if (std::abs(rc[col]) < 1e-9) candidates.push_back(col);
2696 bool refactorize =
false;
2699 for (
int i = 0; i < 10; ++i) {
2702 if (num_pivots >= 5)
break;
2703 if (candidates.empty())
break;
2707 std::uniform_int_distribution<int>(0, candidates.size() - 1)(random_);
2708 const ColIndex entering_col = candidates[index];
2709 std::swap(candidates[index], candidates.back());
2710 candidates.pop_back();
2720 if (reduced_costs_.NeedsBasisRefactorization()) refactorize =
true;
2724 ComputeDirection(entering_col);
2726 RowIndex leaving_row;
2728 bool local_refactorize =
false;
2730 ChooseLeavingVariableRow(entering_col, fake_rc, &local_refactorize,
2731 &leaving_row, &step_length, &target_bound));
2733 if (local_refactorize)
continue;
2735 if (std::abs(step_length) <= 1e-6)
continue;
2736 if (leaving_row !=
kInvalidRow && std::abs(direction_[leaving_row]) < 0.1) {
2739 const Fractional step = (fake_rc > 0.0) ? -step_length : step_length;
2745 const auto get_diff = [
this](ColIndex col,
Fractional old_value,
2747 if (col >= integrality_scale_.size() || integrality_scale_[col] == 0.0) {
2750 const Fractional s = integrality_scale_[col];
2751 return (std::abs(new_value * s - std::round(new_value * s)) -
2752 std::abs(old_value * s - std::round(old_value * s)));
2754 Fractional diff = get_diff(entering_col, variable_values_.Get(entering_col),
2755 variable_values_.Get(entering_col) + step);
2756 for (
const auto e : direction_) {
2757 const ColIndex col = basis_[e.row()];
2758 const Fractional old_value = variable_values_.Get(col);
2759 const Fractional new_value = old_value - e.coefficient() * step;
2760 diff += get_diff(col, old_value, new_value);
2764 if (diff > -1e-2)
continue;
2769 variable_values_.UpdateOnPivoting(direction_, entering_col, step);
2774 SetNonBasicVariableStatusAndDeriveValue(entering_col,
2776 }
else if (step < 0.0) {
2777 SetNonBasicVariableStatusAndDeriveValue(entering_col,
2784 const ColIndex leaving_col = basis_[leaving_row];
2785 update_row_.ComputeUpdateRow(leaving_row);
2794 primal_edge_norms_.UpdateBeforeBasisPivot(
2795 entering_col, leaving_col, leaving_row, direction_, &update_row_);
2796 dual_edge_norms_.UpdateBeforeBasisPivot(
2797 entering_col, leaving_row, direction_,
2798 update_row_.GetUnitRowLeftInverse());
2803 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row, direction_,
2806 const Fractional dir = -direction_[leaving_row] * step;
2807 const bool is_degenerate =
2809 (dir > 0.0 && variable_values_.Get(leaving_col) >= target_bound) ||
2810 (dir < 0.0 && variable_values_.Get(leaving_col) <= target_bound);
2811 if (!is_degenerate) {
2812 variable_values_.Set(leaving_col, target_bound);
2815 UpdateAndPivot(entering_col, leaving_row, target_bound));
2818 VLOG(1) <<
"Polish num_pivots: " << num_pivots <<
" gain:" << total_gain;
2839 Cleanup update_deterministic_time_on_return(
2841 num_consecutive_degenerate_iterations_ = 0;
2842 bool refactorize =
false;
2843 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
2847 primal_prices_.ForceRecomputation();
2849 if (phase_ == Phase::FEASIBILITY) {
2852 objective_.AssignToZero(num_cols_);
2853 variable_values_.UpdatePrimalPhaseICosts(
2854 util::IntegerRange<RowIndex>(RowIndex(0), num_rows_), &objective_);
2855 reduced_costs_.ResetForNewObjective();
2868 if (!refactorize && reduced_costs_.NeedsBasisRefactorization()) {
2869 last_refactorization_reason_ = RefactorizationReason::RC;
2872 if (!refactorize && primal_edge_norms_.NeedsBasisRefactorization()) {
2873 last_refactorization_reason_ = RefactorizationReason::NORM;
2878 if (basis_factorization_.IsRefactorized()) {
2879 CorrectErrorsOnVariableValues();
2880 DisplayIterationInfo(
true, last_refactorization_reason_);
2881 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
2883 if (phase_ == Phase::FEASIBILITY) {
2886 if (variable_values_.UpdatePrimalPhaseICosts(
2887 util::IntegerRange<RowIndex>(RowIndex(0), num_rows_),
2889 reduced_costs_.ResetForNewObjective();
2895 if (phase_ == Phase::OPTIMIZATION &&
2896 ComputeObjectiveValue() < primal_objective_limit_) {
2897 VLOG(1) <<
"Stopping the primal simplex because"
2898 <<
" the objective limit " << primal_objective_limit_
2899 <<
" has been reached.";
2901 objective_limit_reached_ =
true;
2904 }
else if (phase_ == Phase::FEASIBILITY) {
2907 if (variable_values_.UpdatePrimalPhaseICosts(direction_.non_zeros,
2909 reduced_costs_.ResetForNewObjective();
2913 const ColIndex entering_col = primal_prices_.GetBestEnteringColumn();
2915 if (reduced_costs_.AreReducedCostsPrecise() &&
2916 basis_factorization_.IsRefactorized()) {
2917 if (phase_ == Phase::FEASIBILITY) {
2919 variable_values_.ComputeMaximumPrimalInfeasibility();
2920 if (primal_infeasibility <
2921 parameters_.primal_feasibility_tolerance()) {
2924 VLOG(1) <<
"Infeasible problem! infeasibility = "
2925 << primal_infeasibility;
2934 VLOG(1) <<
"Optimal reached, double checking...";
2935 reduced_costs_.MakeReducedCostsPrecise();
2937 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
2941 DCHECK(reduced_costs_.IsValidPrimalEnteringCandidate(entering_col));
2944 ComputeDirection(entering_col);
2949 if (!primal_edge_norms_.TestEnteringEdgeNormPrecision(entering_col,
2951 primal_prices_.RecomputePriceAt(entering_col);
2955 reduced_costs_.TestEnteringReducedCostPrecision(entering_col,
2960 primal_prices_.RecomputePriceAt(entering_col);
2961 if (!reduced_costs_.IsValidPrimalEnteringCandidate(entering_col)) {
2962 reduced_costs_.MakeReducedCostsPrecise();
2963 VLOG(1) <<
"Skipping col #" << entering_col
2964 <<
" whose reduced cost is no longer valid under precise reduced "
2974 if (num_iterations_ == parameters_.max_number_of_iterations())
break;
2977 RowIndex leaving_row;
2979 if (phase_ == Phase::FEASIBILITY) {
2980 PrimalPhaseIChooseLeavingVariableRow(entering_col, reduced_cost,
2981 &refactorize, &leaving_row,
2982 &step_length, &target_bound);
2985 ChooseLeavingVariableRow(entering_col, reduced_cost, &refactorize,
2986 &leaving_row, &step_length, &target_bound));
2989 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
2997 if (!basis_factorization_.IsRefactorized() ||
2998 !reduced_costs_.AreReducedCostsPrecise()) {
2999 VLOG(1) <<
"Infinite step length, double checking...";
3000 reduced_costs_.MakeReducedCostsPrecise();
3002 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3005 if (phase_ == Phase::FEASIBILITY) {
3007 VLOG(1) <<
"Unbounded feasibility problem !?";
3011 solution_primal_ray_.AssignToZero(num_cols_);
3012 for (RowIndex row(0); row < num_rows_; ++row) {
3013 const ColIndex col = basis_[row];
3014 solution_primal_ray_[col] = -direction_[row];
3016 solution_primal_ray_[entering_col] = 1.0;
3017 if (reduced_cost > 0.0) {
3024 Fractional step = (reduced_cost > 0.0) ? -step_length : step_length;
3025 if (phase_ == Phase::FEASIBILITY && leaving_row !=
kInvalidRow) {
3035 step = ComputeStepToMoveBasicVariableToBound(leaving_row, target_bound);
3039 const ColIndex leaving_col =
3045 bool is_degenerate =
false;
3047 Fractional dir = -direction_[leaving_row] * step;
3050 (dir > 0.0 && variable_values_.Get(leaving_col) >= target_bound) ||
3051 (dir < 0.0 && variable_values_.Get(leaving_col) <= target_bound);
3055 if (!is_degenerate) {
3056 DCHECK_EQ(step, ComputeStepToMoveBasicVariableToBound(leaving_row,
3061 variable_values_.UpdateOnPivoting(direction_, entering_col, step);
3064 primal_edge_norms_.UpdateBeforeBasisPivot(
3065 entering_col, basis_[leaving_row], leaving_row, direction_,
3067 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row,
3068 direction_, &update_row_);
3069 primal_prices_.UpdateBeforeBasisPivot(entering_col, &update_row_);
3070 if (!is_degenerate) {
3076 variable_values_.Set(leaving_col, target_bound);
3079 UpdateAndPivot(entering_col, leaving_row, target_bound));
3081 if (is_degenerate) {
3082 timer.AlsoUpdate(&iteration_stats_.degenerate);
3084 timer.AlsoUpdate(&iteration_stats_.normal);
3091 variables_info_.GetTypeRow()[entering_col]);
3093 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3095 }
else if (step < 0.0) {
3096 SetNonBasicVariableStatusAndDeriveValue(entering_col,
3099 primal_prices_.SetAndDebugCheckThatColumnIsDualFeasible(entering_col);
3103 if (phase_ == Phase::FEASIBILITY && leaving_row !=
kInvalidRow) {
3105 variable_values_.SetNonBasicVariableValueFromStatus(leaving_col);
3108 reduced_costs_.SetNonBasicVariableCostToZero(leaving_col,
3109 &objective_[leaving_col]);
3110 primal_prices_.RecomputePriceAt(leaving_col);
3114 if (step_length == 0.0) {
3115 num_consecutive_degenerate_iterations_++;
3117 if (num_consecutive_degenerate_iterations_ > 0) {
3118 iteration_stats_.degenerate_run_size.Add(
3119 num_consecutive_degenerate_iterations_);
3120 num_consecutive_degenerate_iterations_ = 0;
3125 if (num_consecutive_degenerate_iterations_ > 0) {
3126 iteration_stats_.degenerate_run_size.Add(
3127 num_consecutive_degenerate_iterations_);
3143Status RevisedSimplex::DualMinimize(
bool feasibility_phase,
3145 Cleanup update_deterministic_time_on_return(
3147 num_consecutive_degenerate_iterations_ = 0;
3148 bool refactorize =
false;
3149 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
3151 bound_flip_candidates_.clear();
3154 RowIndex leaving_row;
3159 ColIndex entering_col;
3173 const bool old_refactorize_value = refactorize;
3174 if (!refactorize && reduced_costs_.NeedsBasisRefactorization()) {
3175 last_refactorization_reason_ = RefactorizationReason::RC;
3178 if (!refactorize && dual_edge_norms_.NeedsBasisRefactorization()) {
3179 last_refactorization_reason_ = RefactorizationReason::NORM;
3186 if (basis_factorization_.IsRefactorized()) {
3200 if (feasibility_phase || old_refactorize_value) {
3201 reduced_costs_.MakeReducedCostsPrecise();
3212 if (!feasibility_phase) {
3213 MakeBoxedVariableDualFeasible(
3214 variables_info_.GetNonBasicBoxedVariables(),
3216 variable_values_.RecomputeBasicVariableValues();
3217 variable_values_.RecomputeDualPrices(
3218 parameters_.dual_price_prioritize_norm());
3227 if (phase_ == Phase::OPTIMIZATION &&
3229 ComputeObjectiveValue() > dual_objective_limit_) {
3231 "Stopping the dual simplex because"
3232 " the objective limit ",
3233 dual_objective_limit_,
" has been reached.");
3235 objective_limit_reached_ =
true;
3240 DisplayIterationInfo(
false, last_refactorization_reason_);
3241 last_refactorization_reason_ = RefactorizationReason::DEFAULT;
3245 if (!feasibility_phase) {
3248 MakeBoxedVariableDualFeasible(bound_flip_candidates_,
3250 bound_flip_candidates_.clear();
3254 variable_values_.UpdateDualPrices(direction_.non_zeros);
3258 if (feasibility_phase) {
3260 &leaving_row, &cost_variation, &target_bound));
3263 &leaving_row, &cost_variation, &target_bound));
3268 if (!basis_factorization_.IsRefactorized() ||
3269 reduced_costs_.HasCostShift()) {
3270 VLOG(1) <<
"Optimal reached, double checking.";
3271 reduced_costs_.ClearAndRemoveCostShifts();
3274 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3277 if (feasibility_phase) {
3282 if (num_dual_infeasible_positions_ == 0) {
3285 VLOG(1) <<
"DUAL infeasible in dual phase I.";
3295 update_row_.ComputeUnitRowLeftInverse(leaving_row);
3296 if (!dual_edge_norms_.TestPrecision(leaving_row,
3297 update_row_.GetUnitRowLeftInverse())) {
3301 if (feasibility_phase) {
3302 const Fractional price = dual_pricing_vector_[leaving_row];
3304 dual_edge_norms_.GetEdgeSquaredNorms();
3305 dual_prices_.AddOrUpdate(leaving_row,
3306 Square(price) / squared_norms[leaving_row]);
3308 variable_values_.UpdateDualPrices({leaving_row});
3312 update_row_.ComputeUpdateRow(leaving_row);
3314 if (feasibility_phase) {
3316 reduced_costs_.AreReducedCostsPrecise(), update_row_, cost_variation,
3320 reduced_costs_.AreReducedCostsPrecise(), update_row_, cost_variation,
3321 &bound_flip_candidates_, &entering_col));
3326 if (!reduced_costs_.AreReducedCostsPrecise()) {
3327 VLOG(1) <<
"No entering column. Double checking...";
3330 last_refactorization_reason_ = RefactorizationReason::FINAL_CHECK;
3333 DCHECK(basis_factorization_.IsRefactorized());
3334 if (feasibility_phase) {
3336 VLOG(1) <<
"Unbounded dual feasibility problem !?";
3340 solution_dual_ray_ =
3341 Transpose(update_row_.GetUnitRowLeftInverse().values);
3342 update_row_.ComputeFullUpdateRow(leaving_row,
3343 &solution_dual_ray_row_combination_);
3344 if (cost_variation < 0) {
3346 ChangeSign(&solution_dual_ray_row_combination_);
3356 const Fractional entering_coeff = update_row_.GetCoefficient(entering_col);
3357 if (std::abs(entering_coeff) < parameters_.dual_small_pivot_threshold() &&
3358 !reduced_costs_.AreReducedCostsPrecise()) {
3359 VLOG(1) <<
"Trying not to pivot by " << entering_coeff;
3362 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
3366 ComputeDirection(entering_col);
3372 if (std::abs(direction_[leaving_row]) <
3373 parameters_.small_pivot_threshold() * direction_infinity_norm_) {
3374 if (!reduced_costs_.AreReducedCostsPrecise()) {
3375 VLOG(1) <<
"Trying not pivot by " << entering_coeff <<
" ("
3376 << direction_[leaving_row]
3377 <<
") because the direction has a norm of "
3378 << direction_infinity_norm_;
3381 last_refactorization_reason_ = RefactorizationReason::SMALL_PIVOT;
3390 if (std::abs(direction_[leaving_row]) <= 1e-20) {
3391 const std::string error_message = absl::StrCat(
3392 "trying to pivot with number too small: ", direction_[leaving_row]);
3401 if (num_iterations_ == parameters_.max_number_of_iterations()) {
3415 const bool increasing_rc_is_needed =
3416 (cost_variation > 0.0) == (entering_coeff > 0.0);
3417 reduced_costs_.ShiftCostIfNeeded(increasing_rc_is_needed, entering_col);
3420 if (reduced_costs_.StepIsDualDegenerate(increasing_rc_is_needed,
3422 timer.AlsoUpdate(&iteration_stats_.degenerate);
3424 timer.AlsoUpdate(&iteration_stats_.normal);
3433 reduced_costs_.UpdateBeforeBasisPivot(entering_col, leaving_row, direction_,
3435 dual_edge_norms_.UpdateBeforeBasisPivot(
3436 entering_col, leaving_row, direction_,
3437 update_row_.GetUnitRowLeftInverse());
3442 if (feasibility_phase) {
3443 DualPhaseIUpdatePrice(leaving_row, entering_col);
3446 ComputeStepToMoveBasicVariableToBound(leaving_row, target_bound);
3447 variable_values_.UpdateOnPivoting(direction_, entering_col, primal_step);
3451 const ColIndex leaving_col = basis_[leaving_row];
3453 UpdateAndPivot(entering_col, leaving_row, target_bound));
3459 variable_values_.SetNonBasicVariableValueFromStatus(leaving_col);
3468 Cleanup update_deterministic_time_on_return(
3470 bool refactorize =
false;
3474 primal_edge_norms_.Clear();
3475 dual_edge_norms_.Clear();
3476 update_row_.Invalidate();
3477 reduced_costs_.ClearAndRemoveCostShifts();
3479 std::vector<ColIndex> super_basic_cols;
3480 for (
const ColIndex col : variables_info_.GetNotBasicBitRow()) {
3482 variable_values_.Get(col) != 0) {
3483 super_basic_cols.push_back(col);
3487 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 deterministic_random_.seed(parameters.random_seed());
3686 absl_random_ = absl::BitGen(absl::SeedSeq({parameters.random_seed()}));
3687 initial_parameters_ = parameters;
3688 parameters_ = parameters;
3689 PropagateParameters();
3692void RevisedSimplex::PropagateParameters() {
3702void RevisedSimplex::DisplayIterationInfo(
bool primal,
3703 RefactorizationReason reason) {
3705 const std::string first_word = primal ?
"Primal " :
"Dual ";
3712 if (reason != RefactorizationReason::DEFAULT) {
3714 case RefactorizationReason::DEFAULT:
3715 info =
" [default]";
3717 case RefactorizationReason::SMALL_PIVOT:
3718 info =
" [small pivot]";
3720 case RefactorizationReason::IMPRECISE_PIVOT:
3721 info =
" [imprecise pivot]";
3723 case RefactorizationReason::NORM:
3726 case RefactorizationReason::RC:
3727 info =
" [reduced costs]";
3729 case RefactorizationReason::VAR_VALUES:
3730 info =
" [var values]";
3732 case RefactorizationReason::FINAL_CHECK:
3739 case Phase::FEASIBILITY: {
3740 const int64_t iter = num_iterations_;
3743 if (parameters_.use_dual_simplex()) {
3744 if (parameters_.use_dedicated_dual_feasibility_algorithm()) {
3745 objective = reduced_costs_.ComputeSumOfDualInfeasibilities();
3750 objective_,
Transpose(variable_values_.GetDenseRow()));
3752 name =
"sum_dual_infeasibilities";
3754 objective = variable_values_.ComputeSumOfPrimalInfeasibilities();
3755 name =
"sum_primal_infeasibilities";
3758 SOLVER_LOG(logger_, first_word,
"feasibility phase, iteration # ", iter,
3759 ", ", name,
" = ", absl::StrFormat(
"%.15E", objective), info);
3762 case Phase::OPTIMIZATION: {
3763 const int64_t iter = num_iterations_ - num_feasibility_iterations_;
3769 const Fractional objective = ComputeInitialProblemObjectiveValue();
3770 SOLVER_LOG(logger_, first_word,
"optimization phase, iteration # ", iter,
3771 ", objective = ", absl::StrFormat(
"%.15E", objective), info);
3775 const int64_t iter = num_iterations_ - num_feasibility_iterations_ -
3776 num_optimization_iterations_;
3777 SOLVER_LOG(logger_, first_word,
"push phase, iteration # ", iter,
3778 ", remaining_variables_to_push = ",
3779 ComputeNumberOfSuperBasicVariables(), info);
3784void RevisedSimplex::DisplayErrors() {
3785 if (!logger_->LoggingIsEnabled())
return;
3788 SOLVER_LOG(logger_,
"Primal infeasibility (bounds) = ",
3789 variable_values_.ComputeMaximumPrimalInfeasibility());
3790 SOLVER_LOG(logger_,
"Primal residual |A.x - b| = ",
3791 variable_values_.ComputeMaximumPrimalResidual());
3792 SOLVER_LOG(logger_,
"Dual infeasibility (reduced costs) = ",
3793 reduced_costs_.ComputeMaximumDualInfeasibility());
3794 SOLVER_LOG(logger_,
"Dual residual |c_B - y.B| = ",
3795 reduced_costs_.ComputeMaximumDualResidual());
3800std::string StringifyMonomialWithFlags(
const Fractional a,
3801 absl::string_view x) {
3803 a, x, absl::GetFlag(FLAGS_simplex_display_numbers_as_fractions));
3809std::string StringifyWithFlags(
const Fractional x) {
3811 absl::GetFlag(FLAGS_simplex_display_numbers_as_fractions));
3816std::string RevisedSimplex::SimpleVariableInfo(ColIndex col)
const {
3818 VariableType variable_type = variables_info_.GetTypeRow()[col];
3819 VariableStatus variable_status = variables_info_.GetStatusRow()[col];
3820 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
3821 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
3822 absl::StrAppendFormat(&output,
"%d (%s) = %s, %s, %s, [%s,%s]", col.value(),
3823 variable_name_[col],
3824 StringifyWithFlags(variable_values_.Get(col)),
3827 StringifyWithFlags(lower_bounds[col]),
3828 StringifyWithFlags(upper_bounds[col]));
3832void RevisedSimplex::DisplayInfoOnVariables()
const {
3833 if (VLOG_IS_ON(3)) {
3834 for (ColIndex col(0); col < num_cols_; ++col) {
3835 const Fractional variable_value = variable_values_.Get(col);
3836 const Fractional objective_coefficient = objective_[col];
3838 objective_coefficient * variable_value;
3839 VLOG(3) << SimpleVariableInfo(col) <<
". " << variable_name_[col] <<
" = "
3840 << StringifyWithFlags(variable_value) <<
" * "
3841 << StringifyWithFlags(objective_coefficient)
3842 <<
"(obj) = " << StringifyWithFlags(objective_contribution);
3844 VLOG(3) <<
"------";
3848void RevisedSimplex::DisplayVariableBounds() {
3849 if (VLOG_IS_ON(3)) {
3851 const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
3852 const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
3853 for (ColIndex col(0); col < num_cols_; ++col) {
3854 switch (variable_type[col]) {
3858 VLOG(3) << variable_name_[col]
3859 <<
" >= " << StringifyWithFlags(lower_bounds[col]) <<
";";
3862 VLOG(3) << variable_name_[col]
3863 <<
" <= " << StringifyWithFlags(upper_bounds[col]) <<
";";
3866 VLOG(3) << StringifyWithFlags(lower_bounds[col])
3867 <<
" <= " << variable_name_[col]
3868 <<
" <= " << StringifyWithFlags(upper_bounds[col]) <<
";";
3871 VLOG(3) << variable_name_[col] <<
" = "
3872 << StringifyWithFlags(lower_bounds[col]) <<
";";
3875 LOG(DFATAL) <<
"Column " << col <<
" has no meaningful status.";
3882util_intops::StrongVector<RowIndex, SparseRow>
3885 for (ColIndex col(0); col < num_cols_; ++col) {
3886 ComputeDirection(col);
3887 for (
const auto e : direction_) {
3888 if (column_scales ==
nullptr) {
3889 dictionary[e.row()].SetCoefficient(col, e.coefficient());
3893 col < column_scales->
size() ? (*column_scales)[col] : 1.0;
3895 ? (*column_scales)[
GetBasis(e.row())]
3897 dictionary[e.row()].SetCoefficient(
3898 col, direction_[e.row()] * (numerator / denominator));
3907 Status status = Initialize(linear_program);
3909 variable_values_.RecomputeBasicVariableValues();
3910 solution_objective_value_ = ComputeInitialProblemObjectiveValue();
3914void RevisedSimplex::DisplayRevisedSimplexDebugInfo() {
3915 if (VLOG_IS_ON(3)) {
3917 DisplayInfoOnVariables();
3919 std::string output =
"z = " + StringifyWithFlags(ComputeObjectiveValue());
3922 absl::StrAppend(&output, StringifyMonomialWithFlags(reduced_costs[col],
3923 variable_name_[col]));
3925 VLOG(3) << output <<
";";
3927 const RevisedSimplexDictionary dictionary(
nullptr,
this);
3929 for (
const SparseRow& row : dictionary) {
3931 ColIndex basic_col = basis_[r];
3932 absl::StrAppend(&output, variable_name_[basic_col],
" = ",
3933 StringifyWithFlags(variable_values_.
Get(basic_col)));
3934 for (
const SparseRowEntry e : row) {
3935 if (e.col() != basic_col) {
3936 absl::StrAppend(&output,
3937 StringifyMonomialWithFlags(e.coefficient(),
3938 variable_name_[e.col()]));
3941 VLOG(3) << output <<
";";
3943 VLOG(3) <<
"------";
3944 DisplayVariableBounds();
3949void RevisedSimplex::DisplayProblem()
const {
3952 if (VLOG_IS_ON(3)) {
3953 DisplayInfoOnVariables();
3954 std::string output =
"min: ";
3955 bool has_objective =
false;
3956 for (ColIndex col(0); col < num_cols_; ++col) {
3958 has_objective |= (coeff != 0.0);
3959 absl::StrAppend(&output,
3960 StringifyMonomialWithFlags(coeff, variable_name_[col]));
3962 if (!has_objective) {
3963 absl::StrAppend(&output,
" 0");
3965 VLOG(3) << output <<
";";
3966 for (RowIndex row(0); row < num_rows_; ++row) {
3968 for (ColIndex col(0); col < num_cols_; ++col) {
3969 absl::StrAppend(&output,
3970 StringifyMonomialWithFlags(
3971 compact_matrix_.column(col).LookUpCoefficient(row),
3972 variable_name_[col]));
3974 VLOG(3) << output <<
" = 0;";
3976 VLOG(3) <<
"------";
3980void RevisedSimplex::AdvanceDeterministicTime(TimeLimit*
time_limit) {
3983 const double deterministic_time_delta =
3984 current_deterministic_time - last_deterministic_time_update_;
3985 time_limit->AdvanceDeterministicTime(deterministic_time_delta);
3986 last_deterministic_time_update_ = current_deterministic_time;
3989#undef DCHECK_COL_BOUNDS
3990#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.
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)
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)
constexpr double kInfinity
Infinity for type Fractional.
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)
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,...)