From fe66d48cee746408f2093a083b073aa8d60cec4c Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 28 Mar 2017 16:07:49 +0200 Subject: [PATCH] refactor code --- src/glop/basis_representation.cc | 10 +- src/glop/dual_edge_norms.cc | 2 +- src/glop/entering_variable.cc | 16 +- src/glop/initial_basis.cc | 22 +- src/glop/lu_factorization.cc | 6 +- src/glop/markowitz.cc | 8 +- src/glop/preprocessor.cc | 55 +++-- src/glop/preprocessor.h | 5 +- src/glop/primal_edge_norms.cc | 6 +- src/glop/reduced_costs.cc | 14 +- src/glop/revised_simplex.cc | 368 +++++++++++++++---------------- src/glop/revised_simplex.h | 47 +++- src/glop/update_row.cc | 33 ++- src/glop/update_row.h | 4 + src/glop/variables_info.cc | 70 ++++-- src/glop/variables_info.h | 10 +- 16 files changed, 395 insertions(+), 281 deletions(-) diff --git a/src/glop/basis_representation.cc b/src/glop/basis_representation.cc index 24e1080f9e..8b67ba0490 100644 --- a/src/glop/basis_representation.cc +++ b/src/glop/basis_representation.cc @@ -477,6 +477,12 @@ void BasisFactorization::RightSolveForProblemColumn( lu_factorization_.RightSolveLForSparseColumn(matrix_.column(col), d, non_zeros); rank_one_factorization_.RightSolveWithNonZeros(d, non_zeros); + if (col >= right_pool_mapping_.size()) { + // This is needed because when we do an incremental solve with only new + // columns, we still reuse the current factorization without calling + // Refactorize() which would have resized this vector. + right_pool_mapping_.resize(col + 1, kInvalidCol); + } if (non_zeros->empty()) { right_pool_mapping_[col] = right_storage_.AddDenseColumn(*d); } else { @@ -546,7 +552,7 @@ Fractional BasisFactorization::ComputeInverseOneNorm() const { Fractional column_norm = 0.0; // Compute sum_i |inverse_ij|. for (RowIndex row(0); row < num_rows; ++row) { - column_norm += fabs(right_hand_side[row]); + column_norm += std::abs(right_hand_side[row]); } // Compute max_j sum_i |inverse_ij| norm = std::max(norm, column_norm); @@ -566,7 +572,7 @@ Fractional BasisFactorization::ComputeInverseInfinityNorm() const { RightSolve(&right_hand_side); // Compute sum_j |inverse_ij|. for (RowIndex row(0); row < num_rows; ++row) { - row_sum[row] += fabs(right_hand_side[row]); + row_sum[row] += std::abs(right_hand_side[row]); } } // Compute max_i sum_j |inverse_ij| diff --git a/src/glop/dual_edge_norms.cc b/src/glop/dual_edge_norms.cc index ba6d7f0ed6..947691cdf0 100644 --- a/src/glop/dual_edge_norms.cc +++ b/src/glop/dual_edge_norms.cc @@ -61,7 +61,7 @@ void DualEdgeNorms::UpdateBeforeBasisPivot( (sqrt(leaving_squared_norm) - sqrt(old_squared_norm)) / sqrt(leaving_squared_norm); stats_.edge_norms_accuracy.Add(estimated_edge_norms_accuracy); - if (fabs(estimated_edge_norms_accuracy) > + if (std::abs(estimated_edge_norms_accuracy) > parameters_.recompute_edges_norm_threshold()) { VLOG(1) << "Recomputing edge norms: " << sqrt(leaving_squared_norm) << " vs " << sqrt(old_squared_norm); diff --git a/src/glop/entering_variable.cc b/src/glop/entering_variable.cc index 97b8bdcf46..200973e903 100644 --- a/src/glop/entering_variable.cc +++ b/src/glop/entering_variable.cc @@ -201,7 +201,7 @@ Status EnteringVariable::DualChooseEnteringColumn( *entering_col = kInvalidCol; bound_flip_candidates->clear(); Fractional best_coeff = -1.0; - Fractional variation_magnitude = fabs(cost_variation); + Fractional variation_magnitude = std::abs(cost_variation); equivalent_entering_choices_.clear(); while (!breakpoints.empty()) { const ColWithRatio top = breakpoints.front(); @@ -327,7 +327,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn( DCHECK_NE(variable_type[col], VariableType::FIXED_VARIABLE); // Skip if the coeff is too small to be a numerically stable pivot. - if (fabs(update_coefficient[col]) < threshold) continue; + if (std::abs(update_coefficient[col]) < threshold) continue; // We will add ratio * coeff to this column. cost_variation makes sure // the leaving variable will be dual-feasible (its update coeff is @@ -340,7 +340,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn( // Only proceed if there is a transition, note that if reduced_costs[col] // is close to zero, then the variable is supposed to be dual-feasible. - if (fabs(reduced_costs[col]) <= dual_feasibility_tolerance) { + if (std::abs(reduced_costs[col]) <= dual_feasibility_tolerance) { // Continue if the variation goes in the dual-feasible direction. if (coeff > 0 && !can_decrease.IsSet(col)) continue; if (coeff < 0 && !can_increase.IsSet(col)) continue; @@ -357,7 +357,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn( // We are sure there is a transition, add it to the set of breakpoints. breakpoints.push_back( - ColWithRatio(col, fabs(reduced_costs[col]), fabs(coeff))); + ColWithRatio(col, std::abs(reduced_costs[col]), std::abs(coeff))); } // Process the breakpoints in priority order. @@ -372,7 +372,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn( // numerically stable pivot. *entering_col = kInvalidCol; *step = -1.0; - Fractional improvement = fabs(cost_variation); + Fractional improvement = std::abs(cost_variation); while (!breakpoints.empty()) { const ColWithRatio top = breakpoints.front(); @@ -390,7 +390,7 @@ Status EnteringVariable::DualPhaseIChooseEnteringColumn( // improvment, we also render it worse if we keep going in the same // direction. if (can_decrease.IsSet(top.col) && can_increase.IsSet(top.col) && - fabs(reduced_costs[top.col]) > threshold) { + std::abs(reduced_costs[top.col]) > threshold) { improvement -= top.coeff_magnitude; } @@ -444,7 +444,7 @@ void EnteringVariable::DantzigChooseEnteringColumn(ColIndex* entering_col) { *entering_col = kInvalidCol; for (const ColIndex col : reduced_costs_->GetDualInfeasiblePositions()) { if (nested_pricing && !unused_columns_.IsSet(col)) continue; - const Fractional unormalized_price = fabs(reduced_costs[col]); + const Fractional unormalized_price = std::abs(reduced_costs[col]); if (normalize) { if (unormalized_price > best_price * matrix_column_norms[col]) { best_price = unormalized_price / matrix_column_norms[col]; @@ -489,7 +489,7 @@ void EnteringVariable::NormalizedChooseEnteringColumn(ColIndex* entering_col) { *entering_col = col; } } else { - const Fractional positive_reduced_cost = fabs(reduced_costs[col]); + const Fractional positive_reduced_cost = std::abs(reduced_costs[col]); if (positive_reduced_cost >= best_price * weights[col]) { if (positive_reduced_cost == best_price * weights[col]) { equivalent_entering_choices_.push_back(col); diff --git a/src/glop/initial_basis.cc b/src/glop/initial_basis.cc index 42fc7b5965..8e5ac26260 100644 --- a/src/glop/initial_basis.cc +++ b/src/glop/initial_basis.cc @@ -51,7 +51,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, } // This is 'v' in Bixby's paper. - DenseColumn scaled_diagonal_fabs(matrix_.num_rows(), kInfinity); + DenseColumn scaled_diagonal_abs(matrix_.num_rows(), kInfinity); // Compute a list of candidate indices and sort them using the heuristic // described in Bixby's paper. @@ -77,7 +77,7 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, const Fractional kBixbyHighThreshold = 0.99; if (candidate_coeff > kBixbyHighThreshold) { enter_basis = true; - } else if (IsDominated(candidate_col, scaled_diagonal_fabs)) { + } else if (IsDominated(candidate_col, scaled_diagonal_abs)) { candidate_coeff = RestrictedInfinityNorm(candidate_col, can_be_replaced, &candidate_row); if (candidate_coeff != 0.0) { @@ -89,8 +89,8 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols, can_be_replaced[candidate_row] = false; SetSupportToFalse(candidate_col, &has_zero_coefficient); const Fractional kBixbyLowThreshold = 0.01; - scaled_diagonal_fabs[candidate_row] = - kBixbyLowThreshold * fabs(candidate_coeff); + scaled_diagonal_abs[candidate_row] = + kBixbyLowThreshold * std::abs(candidate_coeff); (*basis)[candidate_row] = candidate_col_index; } } @@ -138,7 +138,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols, max_scaled_abs_cost_ = 0.0; for (ColIndex col(0); col < num_cols; ++col) { max_scaled_abs_cost_ = - std::max(max_scaled_abs_cost_, fabs(objective_[col])); + std::max(max_scaled_abs_cost_, std::abs(objective_[col])); if (residual_pattern.ColDegree(col) == 1) { residual_singleton_column.push_back(col); } @@ -172,7 +172,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols, Fractional coeff = 0.0; Fractional max_magnitude = 0.0; for (const SparseColumn::Entry e : matrix_.column(candidate)) { - max_magnitude = std::max(max_magnitude, fabs(e.coefficient())); + max_magnitude = std::max(max_magnitude, std::abs(e.coefficient())); if (can_be_replaced[e.row()]) { row = e.row(); coeff = e.coefficient(); @@ -180,11 +180,11 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols, } } const Fractional kStabilityThreshold = 0.01; - if (fabs(coeff) < kStabilityThreshold * max_magnitude) continue; + if (std::abs(coeff) < kStabilityThreshold * max_magnitude) continue; DCHECK_NE(kInvalidRow, row); partial_diagonal_product *= coeff; - if (fabs(partial_diagonal_product) < kMinimumProductMagnitude) { + if (std::abs(partial_diagonal_product) < kMinimumProductMagnitude) { LOG(INFO) << "Numerical difficulties detected. The product of the " << "diagonal coefficients is currently equal to " << partial_diagonal_product; @@ -204,7 +204,7 @@ bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols, } } - return fabs(partial_diagonal_product) >= kMinimumProductMagnitude; + return std::abs(partial_diagonal_product) >= kMinimumProductMagnitude; } void InitialBasis::ComputeCandidates(ColIndex num_cols, @@ -216,7 +216,7 @@ void InitialBasis::ComputeCandidates(ColIndex num_cols, matrix_.column(col).num_entries() > 0) { candidates->push_back(col); max_scaled_abs_cost_ = - std::max(max_scaled_abs_cost_, fabs(objective_[col])); + std::max(max_scaled_abs_cost_, std::abs(objective_[col])); } } const Fractional kBixbyWeight = 1000.0; @@ -258,7 +258,7 @@ Fractional InitialBasis::GetColumnPenalty(ColIndex col) const { if (type == VariableType::UPPER_AND_LOWER_BOUNDED) { penalty = lower_bound_[col] - upper_bound_[col]; } - return penalty + fabs(objective_[col]) / max_scaled_abs_cost_; + return penalty + std::abs(objective_[col]) / max_scaled_abs_cost_; } bool InitialBasis::BixbyColumnComparator::operator()(ColIndex col_a, diff --git a/src/glop/lu_factorization.cc b/src/glop/lu_factorization.cc index a0f63f4eab..7d35d37d47 100644 --- a/src/glop/lu_factorization.cc +++ b/src/glop/lu_factorization.cc @@ -559,7 +559,7 @@ Fractional LuFactorization::ComputeInverseOneNorm() const { Fractional column_norm = 0.0; // Compute sum_i |basis_matrix_ij|. for (RowIndex row(0); row < num_rows; ++row) { - column_norm += fabs(right_hand_side[row]); + column_norm += std::abs(right_hand_side[row]); } // Compute max_j sum_i |basis_matrix_ij| norm = std::max(norm, column_norm); @@ -579,7 +579,7 @@ Fractional LuFactorization::ComputeInverseInfinityNorm() const { RightSolve(&right_hand_side); // Compute sum_j |basis_matrix_ij|. for (RowIndex row(0); row < num_rows; ++row) { - row_sum[row] += fabs(right_hand_side[row]); + row_sum[row] += std::abs(right_hand_side[row]); } } // Compute max_i sum_j |basis_matrix_ij| @@ -646,7 +646,7 @@ bool LuFactorization::CheckFactorization(const MatrixView& matrix, for (ColIndex col(0); col < should_be_zero.num_cols(); ++col) { for (const SparseColumn::Entry e : should_be_zero.column(col)) { - const Fractional magnitude = fabs(e.coefficient()); + const Fractional magnitude = std::abs(e.coefficient()); if (magnitude > tolerance) { VLOG(2) << magnitude << " != 0, at column " << col; return false; diff --git a/src/glop/markowitz.cc b/src/glop/markowitz.cc index ea01129ae8..0e82ab33ca 100644 --- a/src/glop/markowitz.cc +++ b/src/glop/markowitz.cc @@ -76,7 +76,7 @@ Status Markowitz::ComputeRowAndColumnPermutation(const MatrixView& basis_matrix, // but its magnitude can be really close to zero. In both cases, we // report the singularity of the matrix. if (pivot_row == kInvalidRow || pivot_col == kInvalidCol || - fabs(pivot_coefficient) <= singularity_threshold) { + std::abs(pivot_coefficient) <= singularity_threshold) { RETURN_AND_LOG_ERROR(Status::ERROR_LU, StringPrintf("The matrix is singular! pivot = %E", pivot_coefficient)); @@ -427,7 +427,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm, Fractional max_magnitude = 0.0; for (const SparseColumn::Entry e : column) { - max_magnitude = std::max(max_magnitude, fabs(e.coefficient())); + max_magnitude = std::max(max_magnitude, std::abs(e.coefficient())); } if (max_magnitude == 0.0) { // All symbolic non-zero entries have been cancelled! @@ -438,7 +438,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm, const Fractional skip_threshold = threshold * max_magnitude; for (const SparseColumn::Entry e : column) { - const Fractional magnitude = fabs(e.coefficient()); + const Fractional magnitude = std::abs(e.coefficient()); if (magnitude < skip_threshold) continue; const int row_degree = residual_matrix_non_zero_.RowDegree(e.row()); @@ -446,7 +446,7 @@ int64 Markowitz::FindPivot(const RowPermutation& row_perm, DCHECK_NE(markowitz_number, 0); if (markowitz_number < min_markowitz_number || ((markowitz_number == min_markowitz_number) && - magnitude > fabs(*pivot_coefficient))) { + magnitude > std::abs(*pivot_coefficient))) { min_markowitz_number = markowitz_number; *pivot_col = col; *pivot_row = e.row(); diff --git a/src/glop/preprocessor.cc b/src/glop/preprocessor.cc index e7125ceede..5c5987c84b 100644 --- a/src/glop/preprocessor.cc +++ b/src/glop/preprocessor.cc @@ -306,12 +306,12 @@ VariableStatus ComputeVariableStatus(Fractional value, Fractional lower_bound, // Returns the input with the smallest magnitude or zero if both are infinite. Fractional MinInMagnitudeOrZeroIfInfinite(Fractional a, Fractional b) { - const Fractional value = fabs(a) < fabs(b) ? a : b; + const Fractional value = std::abs(a) < std::abs(b) ? a : b; return IsFinite(value) ? value : 0.0; } Fractional MagnitudeOrZeroIfInfinite(Fractional value) { - return IsFinite(value) ? fabs(value) : 0.0; + return IsFinite(value) ? std::abs(value) : 0.0; } // Returns the maximum magnitude of the finite variable bounds of the given @@ -598,7 +598,7 @@ bool ProportionalColumnPreprocessor::Run(LinearProgram* lp, int num_merged = 0; for (++i; i < sorted_columns.size(); ++i) { if (sorted_columns[i].representative != target_representative) break; - if (fabs(sorted_columns[i].scaled_cost - target_scaled_cost) >= + if (std::abs(sorted_columns[i].scaled_cost - target_scaled_cost) >= parameters_.preprocessor_zero_tolerance()) { break; } @@ -710,7 +710,7 @@ void ProportionalColumnPreprocessor::RecoverSolution( const Fractional bound_factor = column_factors_[col] / column_factors_[representative]; const Fractional scaled_distance = - distance_to_bound[representative] / fabs(bound_factor); + distance_to_bound[representative] / std::abs(bound_factor); const Fractional width = upper_bounds_[col] - lower_bounds_[col]; const bool to_upper_bound = (bound_factor > 0.0) == is_distance_to_upper_bound[representative]; @@ -720,7 +720,7 @@ void ProportionalColumnPreprocessor::RecoverSolution( solution->variable_statuses[col] = ComputeVariableStatus(solution->primal_values[col], lower_bounds_[col], upper_bounds_[col]); - distance_to_bound[representative] -= width * fabs(bound_factor); + distance_to_bound[representative] -= width * std::abs(bound_factor); } else { solution->primal_values[col] = to_upper_bound ? upper_bounds_[col] - scaled_distance @@ -890,7 +890,8 @@ bool ProportionalRowPreprocessor::Run(LinearProgram* lp, // proportionality factor. RowIndex new_representative = lower_source; RowIndex other = upper_source; - if (fabs(row_factors_[new_representative]) < fabs(row_factors_[other])) { + if (std::abs(row_factors_[new_representative]) < + std::abs(row_factors_[other])) { std::swap(new_representative, other); } @@ -1533,7 +1534,7 @@ bool DoubletonFreeColumnPreprocessor::Run(LinearProgram* lp, // We prefer deleting the row with the larger coefficient magnitude because // we will divide by this magnitude. TODO(user): Impact? - if (fabs(r.coeff[DELETED]) < fabs(r.coeff[MODIFIED])) { + if (std::abs(r.coeff[DELETED]) < std::abs(r.coeff[MODIFIED])) { std::swap(r.coeff[DELETED], r.coeff[MODIFIED]); std::swap(r.row[DELETED], r.row[MODIFIED]); } @@ -1577,9 +1578,9 @@ bool DoubletonFreeColumnPreprocessor::Run(LinearProgram* lp, // the numerical error in the formula above, we have a really low // objective instead. The logic is the same as in // AddMultipleToSparseVectorAndIgnoreCommonIndex(). - if (fabs(new_objective) > std::numeric_limits::epsilon() * - 2.0 * - fabs(lp->objective_coefficients()[col])) { + if (std::abs(new_objective) > + std::numeric_limits::epsilon() * 2.0 * + std::abs(lp->objective_coefficients()[col])) { lp->SetObjectiveCoefficient(col, new_objective); } else { lp->SetObjectiveCoefficient(col, 0.0); @@ -1791,10 +1792,10 @@ bool UnconstrainedVariablePreprocessor::Run(LinearProgram* lp, } } - // Change the rhs to reflect the fixed variables. Note that is is important to - // do that after all the calls to RemoveZeroCostUnconstrainedVariable() - // because RemoveZeroCostUnconstrainedVariable() needs to store the rhs before - // this modification! + // Change the rhs to reflect the fixed variables. Note that is important to do + // that after all the calls to RemoveZeroCostUnconstrainedVariable() because + // RemoveZeroCostUnconstrainedVariable() needs to store the rhs before this + // modification! const ColIndex end = column_deletion_helper_.GetMarkedColumns().size(); for (ColIndex col(0); col < end; ++col) { if (column_deletion_helper_.IsColumnMarked(col)) { @@ -1855,7 +1856,7 @@ void UnconstrainedVariablePreprocessor::RecoverSolution( // unbounded direction. if (activity * activity_sign_correction_[row] < 0.0) { const Fractional bound = activity / e.coefficient(); - if (fabs(bound) > fabs(primal_value_shift)) { + if (std::abs(bound) > std::abs(primal_value_shift)) { primal_value_shift = bound; row_at_bound = row; } @@ -2266,7 +2267,7 @@ void SingletonPreprocessor::DeleteSingletonColumnInEquality( // tolerances in a few preprocessors. Like an empty column with a cost of // 1e-17 and unbounded towards infinity is currently implying that the // problem is unbounded. This will need fixing. - if (fabs(new_cost) < parameters_.preprocessor_zero_tolerance()) { + if (std::abs(new_cost) < parameters_.preprocessor_zero_tolerance()) { new_cost = 0.0; } lp->SetObjectiveCoefficient(col, new_cost); @@ -2613,7 +2614,7 @@ bool RemoveNearZeroEntriesPreprocessor::Run(LinearProgram* lp, // and column bounds, and "propagate" the bounds as much as possible so we // can use this better estimate here and remove more near-zero entries. const Fractional max_magnitude = - std::max(fabs(lower_bound), fabs(upper_bound)); + std::max(std::abs(lower_bound), std::abs(upper_bound)); if (max_magnitude == kInfinity || max_magnitude == 0) continue; const Fractional threshold = allowed_impact / max_magnitude; lp->GetMutableSparseColumn(col)->RemoveNearZeroEntriesWithWeights( @@ -2621,7 +2622,7 @@ bool RemoveNearZeroEntriesPreprocessor::Run(LinearProgram* lp, if (lp->objective_coefficients()[col] != 0.0 && num_non_zero_objective_coefficients * - fabs(lp->objective_coefficients()[col]) < + std::abs(lp->objective_coefficients()[col]) < threshold) { lp->SetObjectiveCoefficient(col, 0.0); ++num_zeroed_objective_coefficients; @@ -2756,6 +2757,7 @@ bool DoubletonEqualityRowPreprocessor::Run(LinearProgram* lp, continue; } + // Look at the bounds of both variables and exit early if we can delegate // to another pre-processor; otherwise adjust the bounds of the remaining // variable as necessary. @@ -2777,6 +2779,8 @@ bool DoubletonEqualityRowPreprocessor::Run(LinearProgram* lp, status_ = ProblemStatus::ABNORMAL; break; } + + Fractional lb = r.lb[MODIFIED]; Fractional ub = r.ub[MODIFIED]; Fractional carried_over_lb = @@ -2934,6 +2938,17 @@ void DoubletonEqualityRowPreprocessor::RecoverSolution( } } +void DoubletonEqualityRowPreprocessor:: + SwapDeletedAndModifiedVariableRestoreInfo(RestoreInfo* r) { + using std::swap; + swap(r->col[DELETED], r->col[MODIFIED]); + swap(r->coeff[DELETED], r->coeff[MODIFIED]); + swap(r->lb[DELETED], r->lb[MODIFIED]); + swap(r->ub[DELETED], r->ub[MODIFIED]); + swap(r->column[DELETED], r->column[MODIFIED]); + swap(r->objective_coefficient[DELETED], r->objective_coefficient[MODIFIED]); +} + // -------------------------------------------------------- // DualizerPreprocessor // -------------------------------------------------------- @@ -2949,7 +2964,7 @@ bool DualizerPreprocessor::Run(LinearProgram* lp, TimeLimit* time_limit) { primal_num_rows_ = lp->num_constraints(); primal_is_maximization_problem_ = lp->IsMaximizationProblem(); - // If we need to decide wether or not to take the dual, we only take it when + // If we need to decide whether or not to take the dual, we only take it when // the matrix has more rows than columns. The number of rows of a linear // program gives the size of the square matrices we need to invert and the // order of iterations of the simplex method. So solving a program with less @@ -3219,7 +3234,7 @@ bool ShiftVariableBoundsPreprocessor::Run(LinearProgram* lp, if (0.0 < variable_initial_lbs_[col] || 0.0 > variable_initial_ubs_[col]) { Fractional offset = MinInMagnitudeOrZeroIfInfinite( variable_initial_lbs_[col], variable_initial_ubs_[col]); - if (in_mip_context_ && lp->is_variable_integer()[col]) { + if (in_mip_context_ && lp->IsVariableInteger(col)) { // In the integer case, we truncate the number because if for instance // the lower bound is a positive integer + epsilon, we only want to // shift by the integer and leave the lower bound at epsilon. diff --git a/src/glop/preprocessor.h b/src/glop/preprocessor.h index c3314d7913..9b90dc20ba 100644 --- a/src/glop/preprocessor.h +++ b/src/glop/preprocessor.h @@ -394,7 +394,6 @@ class SingletonPreprocessor : public Preprocessor { ~SingletonPreprocessor() final {} bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final; void RecoverSolution(ProblemSolution* solution) const final; - void UseInMipContext() final { LOG(FATAL) << "Not implemented."; } private: // Returns the MatrixEntry of the given singleton row or column, taking into @@ -550,7 +549,6 @@ class ImpliedFreePreprocessor : public Preprocessor { ~ImpliedFreePreprocessor() final {} bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final; void RecoverSolution(ProblemSolution* solution) const final; - void UseInMipContext() final { LOG(FATAL) << "Not implemented."; } private: // This preprocessor adds fixed offsets to some variables. We remember those @@ -753,7 +751,6 @@ class DoubletonEqualityRowPreprocessor : public Preprocessor { ~DoubletonEqualityRowPreprocessor() final {} bool Run(LinearProgram* linear_program, TimeLimit* time_limit) final; void RecoverSolution(ProblemSolution* solution) const final; - void UseInMipContext() final { LOG(FATAL) << "Not implemented."; } private: enum ColChoice { @@ -799,6 +796,8 @@ class DoubletonEqualityRowPreprocessor : public Preprocessor { ColChoiceAndStatus bound_backtracking_at_lower_bound; ColChoiceAndStatus bound_backtracking_at_upper_bound; }; + void SwapDeletedAndModifiedVariableRestoreInfo(RestoreInfo* r); + std::vector restore_stack_; DISALLOW_COPY_AND_ASSIGN(DoubletonEqualityRowPreprocessor); diff --git a/src/glop/primal_edge_norms.cc b/src/glop/primal_edge_norms.cc index 496e236008..64bb18d91f 100644 --- a/src/glop/primal_edge_norms.cc +++ b/src/glop/primal_edge_norms.cc @@ -80,7 +80,7 @@ void PrimalEdgeNorms::TestEnteringEdgeNormPrecision( const Fractional estimated_edges_norm_accuracy = (precise_norm - sqrt(old_squared_norm)) / precise_norm; stats_.edges_norm_accuracy.Add(estimated_edges_norm_accuracy); - if (fabs(estimated_edges_norm_accuracy) > + if (std::abs(estimated_edges_norm_accuracy) > parameters_.recompute_edges_norm_threshold()) { VLOG(1) << "Recomputing edge norms: " << sqrt(precise_squared_norm) << " vs " << sqrt(old_squared_norm); @@ -285,12 +285,12 @@ void PrimalEdgeNorms::UpdateDevexWeights( // Compared to steepest edge update, the DEVEX weight uses the largest of the // norms of two vectors to approximate the norm of the sum. const Fractional entering_norm = sqrt(PreciseSquaredNorm(direction)); - const Fractional pivot_magnitude = fabs(direction[leaving_row]); + const Fractional pivot_magnitude = std::abs(direction[leaving_row]); const Fractional leaving_norm = std::max(1.0, entering_norm / pivot_magnitude); for (const ColIndex col : update_row.GetNonZeroPositions()) { const Fractional coeff = update_row.GetCoefficient(col); - const Fractional update_vector_norm = fabs(coeff) * leaving_norm; + const Fractional update_vector_norm = std::abs(coeff) * leaving_norm; devex_weights_[col] = std::max(devex_weights_[col], update_vector_norm); } devex_weights_[leaving_col] = leaving_norm; diff --git a/src/glop/reduced_costs.cc b/src/glop/reduced_costs.cc index 39d7550979..d2930ee078 100644 --- a/src/glop/reduced_costs.cc +++ b/src/glop/reduced_costs.cc @@ -95,9 +95,9 @@ bool ReducedCosts::TestEnteringReducedCostPrecision( const Fractional estimated_reduced_costs_accuracy = old_reduced_cost - precise_reduced_cost; const Fractional scale = - (fabs(precise_reduced_cost) <= 1.0) ? 1.0 : precise_reduced_cost; + (std::abs(precise_reduced_cost) <= 1.0) ? 1.0 : precise_reduced_cost; stats_.reduced_costs_accuracy.Add(estimated_reduced_costs_accuracy / scale); - if (fabs(estimated_reduced_costs_accuracy) / scale > + if (std::abs(estimated_reduced_costs_accuracy) / scale > parameters_.recompute_reduced_costs_threshold()) { VLOG(1) << "Recomputing reduced costs: " << precise_reduced_cost << " vs " << old_reduced_cost; @@ -127,7 +127,7 @@ Fractional ReducedCosts::ComputeMaximumDualResidual() const { const Fractional residual = objective_[basic_col] + objective_perturbation_[basic_col] - matrix_.ColumnScalarProduct(basic_col, dual_values); - dual_residual_error = std::max(dual_residual_error, fabs(residual)); + dual_residual_error = std::max(dual_residual_error, std::abs(residual)); } return dual_residual_error; } @@ -144,7 +144,7 @@ Fractional ReducedCosts::ComputeMaximumDualInfeasibility() const { if ((can_increase.IsSet(col) && rc < 0.0) || (can_decrease.IsSet(col) && rc > 0.0)) { maximum_dual_infeasibility = - std::max(maximum_dual_infeasibility, fabs(rc)); + std::max(maximum_dual_infeasibility, std::abs(rc)); } } return maximum_dual_infeasibility; @@ -161,7 +161,7 @@ Fractional ReducedCosts::ComputeSumOfDualInfeasibilities() const { const Fractional rc = reduced_costs_[col]; if ((can_increase.IsSet(col) && rc < 0.0) || (can_decrease.IsSet(col) && rc > 0.0)) { - dual_infeasibility_sum += fabs(fabs(rc)); + dual_infeasibility_sum += std::abs(std::abs(rc)); } } return dual_infeasibility_sum; @@ -372,7 +372,7 @@ void ReducedCosts::ComputeReducedCosts() { // We also compute the dual residual error y.B - c_B. if (is_basic.IsSet(col)) { dual_residual_error = - std::max(dual_residual_error, fabs(reduced_costs_[col])); + std::max(dual_residual_error, std::abs(reduced_costs_[col])); } } } else { @@ -392,7 +392,7 @@ void ReducedCosts::ComputeReducedCosts() { if (is_basic.IsSet(col)) { thread_local_dual_residual_error[omp_get_thread_num()] = std::max(thread_local_dual_residual_error[omp_get_thread_num()], - fabs(reduced_costs_[col])); + std::abs(reduced_costs_[col])); } } // end of omp parallel for diff --git a/src/glop/revised_simplex.cc b/src/glop/revised_simplex.cc index 45de38ac41..02c00a4d5c 100644 --- a/src/glop/revised_simplex.cc +++ b/src/glop/revised_simplex.cc @@ -26,6 +26,7 @@ #include "base/integral_types.h" #include "base/logging.h" #include "base/stringprintf.h" +#include "base/join.h" #include "glop/initial_basis.h" #include "glop/parameters.pb.h" #include "lp_data/lp_data.h" @@ -113,7 +114,7 @@ RevisedSimplex::RevisedSimplex() test_lu_(), feasibility_phase_(true), random_("This is a deterministic seed.") { - PropagateParameters(); + SetParameters(parameters_); } void RevisedSimplex::ClearStateForNextSolve() { @@ -141,7 +142,6 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) { // Initialization. Note That Initialize() must be called first since it // analyzes the current solver state. const double start_time = time_limit->GetElapsedTime(); - random_.Reset(parameters_.random_seed()); RETURN_IF_ERROR(Initialize(lp)); dual_infeasibility_improvement_direction_.clear(); update_row_.Invalidate(); @@ -155,7 +155,7 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) { optimization_time_ = 0.0; total_time_ = 0.0; - if (FLAGS_v > 0) { + if (VLOG_IS_ON(1)) { ComputeNumberOfEmptyRows(); ComputeNumberOfEmptyColumns(); DisplayBasicVariableStatistics(); @@ -174,6 +174,8 @@ Status RevisedSimplex::Solve(const LinearProgram& lp, TimeLimit* time_limit) { current_objective_ = objective_; + // TODO(user): Avoid doing the first phase checks when we know from the + // incremental solve that the solution is already dual or primal feasible. VLOG(1) << "------ First phase: feasibility."; entering_variable_.SetPricingRule(parameters_.feasibility_rule()); if (use_dual) { @@ -421,6 +423,11 @@ const DenseColumn& RevisedSimplex::GetDualRay() const { return solution_dual_ray_; } +const DenseRow& RevisedSimplex::GetDualRayRowCombination() const { + DCHECK_EQ(problem_status_, ProblemStatus::DUAL_UNBOUNDED); + return solution_dual_ray_row_combination_; +} + ColIndex RevisedSimplex::GetBasis(RowIndex row) const { return basis_[row]; } const BasisFactorization& RevisedSimplex::GetBasisFactorization() const { @@ -491,14 +498,14 @@ VariableStatus RevisedSimplex::ComputeDefaultVariableStatus( // Returns the bound with the lowest magnitude. Note that it must be finite // because the VariableStatus::FREE case was tested earlier. DCHECK(IsFinite(lower_bound_[col]) || IsFinite(upper_bound_[col])); - return fabs(lower_bound_[col]) <= fabs(upper_bound_[col]) + return std::abs(lower_bound_[col]) <= std::abs(upper_bound_[col]) ? VariableStatus::AT_LOWER_BOUND : VariableStatus::AT_UPPER_BOUND; } void RevisedSimplex::SetNonBasicVariableStatusAndDeriveValue( ColIndex col, VariableStatus status) { - variables_info_.Update(col, status); + variables_info_.UpdateToNonBasicStatus(col, status); variable_values_.SetNonBasicVariableValueFromStatus(col); } @@ -594,9 +601,9 @@ void RevisedSimplex::UseSingletonColumnInInitialBasis(RowToColMapping* basis) { const Fractional slope = matrix_with_slack_.column(col).GetFirstCoefficient(); if (variable_values_.Get(col) == lower_bound_[col]) { - cost_variation[col] = objective_[col] / fabs(slope); + cost_variation[col] = objective_[col] / std::abs(slope); } else { - cost_variation[col] = -objective_[col] / fabs(slope); + cost_variation[col] = -objective_[col] / std::abs(slope); } singleton_column.push_back(col); } @@ -882,40 +889,31 @@ void RevisedSimplex::InitializeVariableStatusesForWarmStart( status = state.statuses[col - num_new_cols]; } - // Remove incompatibilities between the warm status and the variable bounds. - // We use the default status as an indication of the bounds type. - if ((status == VariableStatus::FREE && - default_status != VariableStatus::FREE) || - (status == VariableStatus::FIXED_VALUE && - default_status != VariableStatus::FIXED_VALUE) || - (status == VariableStatus::AT_LOWER_BOUND && - lower_bound_[col] == -kInfinity) || - (status == VariableStatus::AT_UPPER_BOUND && - upper_bound_[col] == kInfinity)) { - status = default_status; - } - if (default_status == VariableStatus::FIXED_VALUE && - (status == VariableStatus::AT_UPPER_BOUND || - status == VariableStatus::AT_LOWER_BOUND)) { - status = VariableStatus::FIXED_VALUE; - } - - // Do not allow more than num_rows_ VariableStatus::BASIC variables. if (status == VariableStatus::BASIC) { + // Do not allow more than num_rows_ VariableStatus::BASIC variables. if (num_basic_variables == num_rows_) { VLOG(1) << "Too many basic variables in the warm-start basis." << "Only keeping the first ones as VariableStatus::BASIC."; - status = default_status; + SetNonBasicVariableStatusAndDeriveValue(col, default_status); } else { ++num_basic_variables; + variables_info_.UpdateToBasicStatus(col); } - } - - // Set the status. - if (status != VariableStatus::BASIC) { - SetNonBasicVariableStatusAndDeriveValue(col, status); } else { - variables_info_.Update(col, VariableStatus::BASIC); + // Remove incompatibilities between the warm status and the variable + // bounds. We use the default status as an indication of the bounds + // type. + if ((status != default_status) && + ((default_status == VariableStatus::FIXED_VALUE) || + (status == VariableStatus::FREE) || + (status == VariableStatus::FIXED_VALUE) || + (status == VariableStatus::AT_LOWER_BOUND && + lower_bound_[col] == -kInfinity) || + (status == VariableStatus::AT_UPPER_BOUND && + upper_bound_[col] == kInfinity))) { + status = default_status; + } + SetNonBasicVariableStatusAndDeriveValue(col, status); } } } @@ -1083,7 +1081,7 @@ Status RevisedSimplex::Initialize(const LinearProgram& lp) { // Computes the variable name as soon as possible for logging. // TODO(user): do we really need to store them? we could just compute them // on the fly since we do not need the speed. - if (FLAGS_v > 0) { + if (VLOG_IS_ON(1)) { SetVariableNames(); } @@ -1329,7 +1327,7 @@ void RevisedSimplex::ComputeDirection(ColIndex col) { direction_infinity_norm_ = 0.0; for (const RowIndex row : direction_non_zero_) { direction_infinity_norm_ = - std::max(direction_infinity_norm_, fabs(direction_[row])); + std::max(direction_infinity_norm_, std::abs(direction_[row])); } IF_STATS_ENABLED(ratio_test_stats_.direction_density.Add( num_rows_ == 0 ? 0.0 @@ -1352,7 +1350,7 @@ void RevisedSimplex::SkipVariableForRatioTest(RowIndex row) { // during the ratio test. The direction vector will be restored later from // the information in direction_ignored_position_. IF_STATS_ENABLED( - ratio_test_stats_.abs_skipped_pivot.Add(fabs(direction_[row]))); + ratio_test_stats_.abs_skipped_pivot.Add(std::abs(direction_[row]))); VLOG(1) << "Skipping leaving variable with coefficient " << direction_[row]; direction_ignored_position_.SetCoefficient(row, direction_[row]); direction_[row] = 0.0; @@ -1397,7 +1395,7 @@ Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates( leaving_candidates->Clear(); const Fractional threshold = parameters_.ratio_test_zero_threshold(); for (const RowIndex row : direction_non_zero_) { - const Fractional magnitude = fabs(direction_[row]); + const Fractional magnitude = std::abs(direction_[row]); if (magnitude < threshold) continue; Fractional ratio = GetRatio(row); // TODO(user): The perturbation is currently disabled, so no need to test @@ -1405,7 +1403,7 @@ Fractional RevisedSimplex::ComputeHarrisRatioAndLeavingCandidates( if (false && ratio < 0.0) { // If the variable is already pass its bound, we use the perturbed version // of the bound (if bound_perturbation_[basis_[row]] is not zero). - ratio += fabs(bound_perturbation_[basis_[row]] / direction_[row]); + ratio += std::abs(bound_perturbation_[basis_[row]] / direction_[row]); } if (ratio <= harris_ratio) { leaving_candidates->SetCoefficient(row, ratio); @@ -1510,7 +1508,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow( // If the magnitudes are the same, we choose the leaving variable with // what is probably the more stable ratio, see // IsRatioMoreOrEquallyStable(). - const Fractional candidate_magnitude = fabs(direction_[row]); + const Fractional candidate_magnitude = std::abs(direction_[row]); if (candidate_magnitude < pivot_magnitude) continue; if (candidate_magnitude == pivot_magnitude) { if (!IsRatioMoreOrEquallyStable(ratio, current_ratio)) continue; @@ -1590,7 +1588,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow( // Note(user): This reduces quite a bit the number of iterations. // What is its impact? Is it dangerous? - if (fabs(direction_[*leaving_row]) < + if (std::abs(direction_[*leaving_row]) < parameters_.minimum_acceptable_pivot()) { SkipVariableForRatioTest(*leaving_row); continue; @@ -1633,7 +1631,7 @@ Status RevisedSimplex::ChooseLeavingVariableRow( equivalent_leaving_choices_.size()); } if (*leaving_row != kInvalidRow) { - ratio_test_stats_.abs_used_pivot.Add(fabs(direction_[*leaving_row])); + ratio_test_stats_.abs_used_pivot.Add(std::abs(direction_[*leaving_row])); } }); return Status::OK; @@ -1725,7 +1723,7 @@ void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow( for (RowIndex row : direction_non_zero_) { const Fractional direction = reduced_cost > 0.0 ? direction_[row] : -direction_[row]; - const Fractional magnitude = fabs(direction); + const Fractional magnitude = std::abs(direction); if (magnitude < tolerance) continue; // Computes by how much we can add 'direction' to the basic variable value @@ -1766,7 +1764,7 @@ void RevisedSimplex::PrimalPhaseIChooseLeavingVariableRow( // Select the last breakpoint that still improves the infeasibility and has // the largest coefficient magnitude. - Fractional improvement = fabs(reduced_cost); + Fractional improvement = std::abs(reduced_cost); Fractional best_magnitude = 0.0; *leaving_row = kInvalidRow; while (!breakpoints.empty()) { @@ -2635,15 +2633,22 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) { } else { problem_status_ = ProblemStatus::DUAL_UNBOUNDED; solution_dual_ray_ = update_row_.GetUnitRowLeftInverse().dense_column; + update_row_.RecomputeFullUpdateRow(leaving_row); + solution_dual_ray_row_combination_.assign(num_cols_, 0.0); + for (const ColIndex col : update_row_.GetNonZeroPositions()) { + solution_dual_ray_row_combination_[col] = + update_row_.GetCoefficient(col); + } if (cost_variation < 0) { ChangeSign(&solution_dual_ray_); + ChangeSign(&solution_dual_ray_row_combination_); } } return Status::OK; } // If the coefficient is too small, we recompute the reduced costs. - if (fabs(entering_coeff) < parameters_.dual_small_pivot_threshold() && + if (std::abs(entering_coeff) < parameters_.dual_small_pivot_threshold() && !reduced_costs_.AreReducedCostsPrecise()) { VLOG(1) << "Trying not to pivot by " << entering_coeff; refactorize = true; @@ -2655,7 +2660,7 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) { // direction_[leaving_row] is actually 0 which causes a floating // point exception below. ComputeDirection(entering_col); - if (fabs(direction_[leaving_row]) < + if (std::abs(direction_[leaving_row]) < parameters_.minimum_acceptable_pivot()) { VLOG(1) << "Do not pivot by " << entering_coeff << " because the direction is " << direction_[leaving_row]; @@ -2719,7 +2724,8 @@ Status RevisedSimplex::DualMinimize(TimeLimit* time_limit) { // This is slow, but otherwise we have a really bad precision on the // variable values ... - if (fabs(primal_step) * parameters_.primal_feasibility_tolerance() > 1.0) { + if (std::abs(primal_step) * parameters_.primal_feasibility_tolerance() > + 1.0) { refactorize = true; } ++num_iterations_; @@ -2770,6 +2776,7 @@ Fractional RevisedSimplex::ComputeInitialProblemObjectiveValue() const { void RevisedSimplex::SetParameters(const GlopParameters& parameters) { SCOPED_TIME_STAT(&function_stats_); + random_.Reset(parameters_.random_seed()); initial_parameters_ = parameters; parameters_ = parameters; PropagateParameters(); @@ -2788,35 +2795,37 @@ void RevisedSimplex::PropagateParameters() { } void RevisedSimplex::DisplayIterationInfo() const { - if (FLAGS_v < 1) return; - const int iter = feasibility_phase_ - ? num_iterations_ - : num_iterations_ - num_feasibility_iterations_; - // Note that in the dual phase II, ComputeObjectiveValue() is also computing - // the dual objective even if it uses the variable values. This is because if - // we modify the bounds to make the problem primal-feasible, we are at the - // optimal and hence the two objectives are the same. - const Fractional objective = - !feasibility_phase_ - ? ComputeInitialProblemObjectiveValue() - : (parameters_.use_dual_simplex() - ? reduced_costs_.ComputeSumOfDualInfeasibilities() - : variable_values_.ComputeSumOfPrimalInfeasibilities()); - VLOG(1) << (feasibility_phase_ ? "Feasibility" : "Optimization") - << " phase, iteration # " << iter - << ", objective = " << StringPrintf("%.15E", objective); + if (VLOG_IS_ON(1)) { + const int iter = feasibility_phase_ + ? num_iterations_ + : num_iterations_ - num_feasibility_iterations_; + // Note that in the dual phase II, ComputeObjectiveValue() is also computing + // the dual objective even if it uses the variable values. This is because + // if we modify the bounds to make the problem primal-feasible, we are at + // the optimal and hence the two objectives are the same. + const Fractional objective = + !feasibility_phase_ + ? ComputeInitialProblemObjectiveValue() + : (parameters_.use_dual_simplex() + ? reduced_costs_.ComputeSumOfDualInfeasibilities() + : variable_values_.ComputeSumOfPrimalInfeasibilities()); + VLOG(1) << (feasibility_phase_ ? "Feasibility" : "Optimization") + << " phase, iteration # " << iter + << ", objective = " << StringPrintf("%.15E", objective); + } } void RevisedSimplex::DisplayErrors() const { - if (FLAGS_v < 1) return; - VLOG(1) << "Primal infeasibility (bounds) = " - << variable_values_.ComputeMaximumPrimalInfeasibility(); - VLOG(1) << "Primal residual |A.x - b| = " - << variable_values_.ComputeMaximumPrimalResidual(); - VLOG(1) << "Dual infeasibility (reduced costs) = " - << reduced_costs_.ComputeMaximumDualInfeasibility(); - VLOG(1) << "Dual residual |c_B - y.B| = " - << reduced_costs_.ComputeMaximumDualResidual(); + if (VLOG_IS_ON(1)) { + VLOG(1) << "Primal infeasibility (bounds) = " + << variable_values_.ComputeMaximumPrimalInfeasibility(); + VLOG(1) << "Primal residual |A.x - b| = " + << variable_values_.ComputeMaximumPrimalResidual(); + VLOG(1) << "Dual infeasibility (reduced costs) = " + << reduced_costs_.ComputeMaximumDualInfeasibility(); + VLOG(1) << "Dual residual |c_B - y.B| = " + << reduced_costs_.ComputeMaximumDualResidual(); + } } namespace { @@ -2848,144 +2857,129 @@ std::string RevisedSimplex::SimpleVariableInfo(ColIndex col) const { } void RevisedSimplex::DisplayInfoOnVariables() const { - if (FLAGS_v < 3) return; - for (ColIndex col(0); col < num_cols_; ++col) { - const Fractional variable_value = variable_values_.Get(col); - const Fractional objective_coefficient = current_objective_[col]; - const Fractional objective_contribution = - objective_coefficient * variable_value; - VLOG(3) << SimpleVariableInfo(col) << ". " << variable_name_[col] << " = " - << StringifyWithFlags(variable_value) << " * " - << StringifyWithFlags(objective_coefficient) - << "(obj) = " << StringifyWithFlags(objective_contribution); + if (VLOG_IS_ON(3)) { + for (ColIndex col(0); col < num_cols_; ++col) { + const Fractional variable_value = variable_values_.Get(col); + const Fractional objective_coefficient = current_objective_[col]; + const Fractional objective_contribution = + objective_coefficient * variable_value; + VLOG(3) << SimpleVariableInfo(col) << ". " << variable_name_[col] << " = " + << StringifyWithFlags(variable_value) << " * " + << StringifyWithFlags(objective_coefficient) + << "(obj) = " << StringifyWithFlags(objective_contribution); + } + VLOG(3) << "------"; } - VLOG(3) << "------"; } void RevisedSimplex::DisplayVariableBounds() { - if (FLAGS_v < 3) return; - const VariableTypeRow& variable_type = variables_info_.GetTypeRow(); - for (ColIndex col(0); col < num_cols_; ++col) { - switch (variable_type[col]) { - case VariableType::UNCONSTRAINED: - break; - case VariableType::LOWER_BOUNDED: - VLOG(3) << variable_name_[col] - << " >= " << StringifyWithFlags(lower_bound_[col]) << ";"; - break; - case VariableType::UPPER_BOUNDED: - VLOG(3) << variable_name_[col] - << " <= " << StringifyWithFlags(upper_bound_[col]) << ";"; - break; - case VariableType::UPPER_AND_LOWER_BOUNDED: - VLOG(3) << StringifyWithFlags(lower_bound_[col]) - << " <= " << variable_name_[col] - << " <= " << StringifyWithFlags(upper_bound_[col]) << ";"; - break; - case VariableType::FIXED_VARIABLE: - VLOG(3) << variable_name_[col] << " = " - << StringifyWithFlags(lower_bound_[col]) << ";"; - break; - default: // This should never happen. - LOG(DFATAL) << "Column " << col - << " has no meaningful status."; - break; + if (VLOG_IS_ON(3)) { + const VariableTypeRow& variable_type = variables_info_.GetTypeRow(); + for (ColIndex col(0); col < num_cols_; ++col) { + switch (variable_type[col]) { + case VariableType::UNCONSTRAINED: + break; + case VariableType::LOWER_BOUNDED: + VLOG(3) << variable_name_[col] + << " >= " << StringifyWithFlags(lower_bound_[col]) << ";"; + break; + case VariableType::UPPER_BOUNDED: + VLOG(3) << variable_name_[col] + << " <= " << StringifyWithFlags(upper_bound_[col]) << ";"; + break; + case VariableType::UPPER_AND_LOWER_BOUNDED: + VLOG(3) << StringifyWithFlags(lower_bound_[col]) + << " <= " << variable_name_[col] + << " <= " << StringifyWithFlags(upper_bound_[col]) << ";"; + break; + case VariableType::FIXED_VARIABLE: + VLOG(3) << variable_name_[col] << " = " + << StringifyWithFlags(lower_bound_[col]) << ";"; + break; + default: // This should never happen. + LOG(DFATAL) << "Column " << col + << " has no meaningful status."; + break; + } } } } -namespace { -struct MatrixElementForDisplayDictionary { - MatrixElementForDisplayDictionary(ColIndex _col, RowIndex _row, - Fractional _value) - : col(_col), row(_row), value(_value) {} - ColIndex col; - RowIndex row; - Fractional value; -}; - -Fractional DictionaryLookUpValue( - const std::vector& matrix, ColIndex col, - RowIndex row) { - const size_t num_elements = matrix.size(); - for (int i = 0; i < num_elements; ++i) { - if (matrix[i].col == col && matrix[i].row == row) { - return matrix[i].value; - } - } - return Fractional(0.0); -} -} // anonymous namespace - -void RevisedSimplex::DisplayDictionary() { - // This function has a complexity in - // O(num_cols_ * num_rows_ * non zeros in column). - // It uses LookUpValue(row, col) which has a time complexity - // in O(non zeros in column). - // This choice was made because it is not expected that people use - // this function on large problems and/or very dense problems. - // TODO(user): re-evaluate this design decision if need be. - if (FLAGS_v < 4) return; - std::vector dictionary; - - DisplayInfoOnVariables(); +ITIVector RevisedSimplex::ComputeDictionary() { + ITIVector dictionary(num_rows_.value()); for (ColIndex col(0); col < num_cols_; ++col) { ComputeDirection(col); for (const RowIndex row : direction_non_zero_) { - dictionary.push_back( - MatrixElementForDisplayDictionary(col, row, direction_[row])); + dictionary[row].SetCoefficient(col, direction_[row]); } } + return dictionary; +} - std::string output = "z = " + StringifyWithFlags(ComputeObjectiveValue()); - const DenseRow& reduced_costs = reduced_costs_.GetReducedCosts(); - for (const ColIndex col : variables_info_.GetNotBasicBitRow()) { - output += - StringifyMonomialWithFlags(reduced_costs[col], variable_name_[col]); - } - VLOG(3) << output << ";"; - for (RowIndex row(0); row < num_rows_; ++row) { - output = ""; - ColIndex basic_col = basis_[row]; - output += variable_name_[basic_col] + " = " + - StringifyWithFlags(variable_values_.Get(basic_col)); - for (ColIndex col(0); col < num_cols_; ++col) { - if (col != basic_col) { - const Fractional value = DictionaryLookUpValue(dictionary, col, row); - output += StringifyMonomialWithFlags(value, variable_name_[col]); - } +void RevisedSimplex::DisplayRevisedSimplexDebugInfo() { + if (VLOG_IS_ON(3)) { + // This function has a complexity in O(num_non_zeros_in_matrix). + DisplayInfoOnVariables(); + + std::string output = "z = " + StringifyWithFlags(ComputeObjectiveValue()); + const DenseRow& reduced_costs = reduced_costs_.GetReducedCosts(); + for (const ColIndex col : variables_info_.GetNotBasicBitRow()) { + StrAppend(&output, StringifyMonomialWithFlags(reduced_costs[col], + variable_name_[col])); } VLOG(3) << output << ";"; + + const RevisedSimplexDictionary dictionary(this); + RowIndex r(0); + for (const SparseRow& row : dictionary) { + output.clear(); + ColIndex basic_col = basis_[r]; + StrAppend(&output, variable_name_[basic_col], " = ", + StringifyWithFlags(variable_values_.Get(basic_col))); + for (const SparseRowEntry e : row) { + if (e.col() != basic_col) { + StrAppend(&output, + StringifyMonomialWithFlags(e.coefficient(), + variable_name_[e.col()])); + } + } + VLOG(3) << output << ";"; + } + VLOG(3) << "------"; + DisplayVariableBounds(); + ++r; } - VLOG(3) << "------"; - DisplayVariableBounds(); } void RevisedSimplex::DisplayProblem() const { - // TODO(user): see comment for DisplayDictionary() - if (FLAGS_v < 3) return; - DisplayInfoOnVariables(); - std::string output = "min: "; - bool has_objective = false; - for (ColIndex col(0); col < num_cols_; ++col) { - const Fractional coeff = current_objective_[col]; - has_objective |= (coeff != 0.0); - output += StringifyMonomialWithFlags(coeff, variable_name_[col]); - } - if (!has_objective) { - output += " 0"; - } - VLOG(3) << output << ";"; - for (RowIndex row(0); row < num_rows_; ++row) { - output = ""; + // This function has a complexity in O(num_rows * num_cols * + // num_non_zeros_in_row). + if (VLOG_IS_ON(3)) { + DisplayInfoOnVariables(); + std::string output = "min: "; + bool has_objective = false; for (ColIndex col(0); col < num_cols_; ++col) { - output += StringifyMonomialWithFlags( - matrix_with_slack_.column(col).LookUpCoefficient(row), - variable_name_[col]); + const Fractional coeff = current_objective_[col]; + has_objective |= (coeff != 0.0); + StrAppend(&output, + StringifyMonomialWithFlags(coeff, variable_name_[col])); } - VLOG(3) << output << " = 0;"; + if (!has_objective) { + StrAppend(&output, " 0"); + } + VLOG(3) << output << ";"; + for (RowIndex row(0); row < num_rows_; ++row) { + output = ""; + for (ColIndex col(0); col < num_cols_; ++col) { + StrAppend( + &output, StringifyMonomialWithFlags( + matrix_with_slack_.column(col).LookUpCoefficient(row), + variable_name_[col])); + } + VLOG(3) << output << " = 0;"; + } + VLOG(3) << "------"; } - VLOG(3) << "------"; } void RevisedSimplex::AdvanceDeterministicTime(TimeLimit* time_limit) { diff --git a/src/glop/revised_simplex.h b/src/glop/revised_simplex.h index d043a7b6b5..a590802a2d 100644 --- a/src/glop/revised_simplex.h +++ b/src/glop/revised_simplex.h @@ -114,6 +114,7 @@ #include "lp_data/lp_print_utils.h" #include "lp_data/lp_types.h" #include "lp_data/matrix_scaler.h" +#include "lp_data/sparse_row.h" #include "base/random.h" #include "util/time_limit.h" @@ -206,6 +207,9 @@ class RevisedSimplex { const DenseRow& GetPrimalRay() const; const DenseColumn& GetDualRay() const; + // This is the "dual ray" linear combination of the matrix rows. + const DenseRow& GetDualRayRowCombination() const; + // Returns the index of the column in the basis and the basis factorization. // Note that the order of the column in the basis is important since it is the // one used by the various solve functions provided by the BasisFactorization @@ -216,6 +220,11 @@ class RevisedSimplex { // Returns statistics about this class as a std::string. std::string StatString(); + // Computes the dictionary B^-1*N on-the-fly row by row. Returns the resulting + // matrix as a vector of sparse rows so that it is easy to use it on the left + // side in the matrix multiplication. Runs in O(num_non_zeros_in_matrix). + RowMajorSparseMatrix ComputeDictionary(); + private: // Propagates parameters_ to all the other classes that need it. // @@ -254,8 +263,12 @@ class RevisedSimplex { // Displays the bounds of the variables. void DisplayVariableBounds(); - // Displays the Linear Programming problem as a dictionary, - // taking into account the iterations that have been made. + // Displays the following information: + // * Linear Programming problem as a dictionary, taking into + // account the iterations that have been made; + // * Variable info; + // * Reduced costs; + // * Variable bounds. // A dictionary is in the form: // xB = value + sum_{j in N} pa_ij x_j // z = objective_value + sum_{i in N} rc_i x_i @@ -264,7 +277,7 @@ class RevisedSimplex { // after the pivotings. // Dictionaries are the modern way of presenting the result of an iteration // of the Simplex algorithm in the literature. - void DisplayDictionary(); + void DisplayRevisedSimplexDebugInfo(); // Displays the Linear Programming problem as it was input. void DisplayProblem() const; @@ -385,7 +398,7 @@ class RevisedSimplex { // See Chvatal's book for more detail (Chapter 7). void ComputeDirection(ColIndex col); - // Computes a - B.d in error_ and return the maximum fabs() of its coeffs. + // Computes a - B.d in error_ and return the maximum std::abs() of its coeffs. // TODO(user): Use this to trigger a refactorization of B? Or to track the // error created by calls to SkipVariableForRatioTest()? Fractional ComputeDirectionError(ColIndex col); @@ -633,6 +646,7 @@ class RevisedSimplex { DenseRow solution_reduced_costs_; DenseRow solution_primal_ray_; DenseColumn solution_dual_ray_; + DenseRow solution_dual_ray_row_combination_; BasisState solution_state_; bool solution_state_has_been_set_externally_; @@ -771,6 +785,31 @@ class RevisedSimplex { DISALLOW_COPY_AND_ASSIGN(RevisedSimplex); }; +// Hides the details of the dictionary matrix implementation. In the future, +// GLOP will support generating the dictionary one row at a time without having +// to store the whole matrix in memory. +class RevisedSimplexDictionary { + public: + typedef RowMajorSparseMatrix::const_iterator ConstIterator; + + // RevisedSimplex cannot be passed const because we have to call a non-const + // method ComputeDictionary. + explicit RevisedSimplexDictionary(RevisedSimplex* revised_simplex) + : dictionary_(CHECK_NOTNULL(revised_simplex)->ComputeDictionary()) {} + + ConstIterator begin() const { return dictionary_.begin(); } + ConstIterator end() const { return dictionary_.end(); } + + size_t size() const { return dictionary_.size(); } + + private: + const RowMajorSparseMatrix dictionary_; + DISALLOW_COPY_AND_ASSIGN(RevisedSimplexDictionary); +}; + +// TODO(user): When a row-by-row generation of the dictionary is supported, +// implement DictionaryIterator class that would call it inside operator*(). + } // namespace glop } // namespace operations_research diff --git a/src/glop/update_row.cc b/src/glop/update_row.cc index 8419f9dd04..3c8266e516 100644 --- a/src/glop/update_row.cc +++ b/src/glop/update_row.cc @@ -162,7 +162,7 @@ void UpdateRow::ComputeUpdatesRowWise() { non_zero_position_list_.clear(); const Fractional drop_tolerance = parameters_.drop_tolerance(); for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) { - if (fabs(coefficient_[col]) > drop_tolerance) { + if (std::abs(coefficient_[col]) > drop_tolerance) { non_zero_position_list_.push_back(col); } } @@ -202,12 +202,37 @@ void UpdateRow::ComputeUpdatesRowWiseHypersparse() { // non-zero coefficients contiguously in a vector is better than keeping // them as they are. Note however that we will iterate only twice on the // update row coefficients during an iteration. - if (fabs(coefficient_[col]) > drop_tolerance) { + if (std::abs(coefficient_[col]) > drop_tolerance) { non_zero_position_list_.push_back(col); } } } +// Note that we use the same algo as ComputeUpdatesColumnWise() here. The +// others version might be faster, but this is called only once per solve, so +// it shouldn't be too bad. +void UpdateRow::RecomputeFullUpdateRow(RowIndex leaving_row) { + CHECK(!compute_update_row_); + const ColIndex num_cols = matrix_.num_cols(); + const Fractional drop_tolerance = parameters_.drop_tolerance(); + coefficient_.resize(num_cols, 0.0); + non_zero_position_list_.clear(); + + // Fills the only position at one in the basic columns. + coefficient_[basis_[leaving_row]] = 1.0; + non_zero_position_list_.push_back(basis_[leaving_row]); + + // Fills the non-basic column. + for (const ColIndex col : variables_info_.GetNotBasicBitRow()) { + const Fractional coeff = + matrix_.ColumnScalarProduct(col, unit_row_left_inverse_); + if (std::abs(coeff) > drop_tolerance) { + non_zero_position_list_.push_back(col); + coefficient_[col] = coeff; + } + } +} + void UpdateRow::ComputeUpdatesColumnWise() { SCOPED_TIME_STAT(&stats_); @@ -227,7 +252,7 @@ void UpdateRow::ComputeUpdatesColumnWise() { // sparsity. Note that it shouldn't be too bad to use a non-zero drop // tolerance here because even if we introduce some precision issues, the // quantities updated by this update row will eventually be recomputed. - if (fabs(coeff) > drop_tolerance) { + if (std::abs(coeff) > drop_tolerance) { non_zero_position_list_.push_back(col); coefficient_[col] = coeff; } @@ -248,7 +273,7 @@ void UpdateRow::ComputeUpdatesColumnWise() { } // End of omp parallel for. for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) { - if (fabs(coefficient_[col]) > drop_tolerance) { + if (std::abs(coefficient_[col]) > drop_tolerance) { non_zero_position_list_.push_back(col); } } diff --git a/src/glop/update_row.h b/src/glop/update_row.h index 9928f3c104..f7b3b49a26 100644 --- a/src/glop/update_row.h +++ b/src/glop/update_row.h @@ -66,6 +66,10 @@ class UpdateRow { return coefficient_[col]; } + // This must be called after a call to ComputeUpdateRow(). It will fill + // all the non-relevant positions that where not filled by ComputeUpdateRow(). + void RecomputeFullUpdateRow(RowIndex leaving_row); + // Sets to zero the coefficient for column col. void IgnoreUpdatePosition(ColIndex col); diff --git a/src/glop/variables_info.cc b/src/glop/variables_info.cc index 5dc7455cd6..26caabb950 100644 --- a/src/glop/variables_info.cc +++ b/src/glop/variables_info.cc @@ -32,7 +32,7 @@ void VariablesInfo::Initialize() { variable_status_.resize(num_cols, VariableStatus::FREE); can_increase_.ClearAndResize(num_cols); can_decrease_.ClearAndResize(num_cols); - is_relevant_.ClearAndResize(num_cols); + relevance_.ClearAndResize(num_cols); is_basic_.ClearAndResize(num_cols); not_basic_.ClearAndResize(num_cols); non_basic_boxed_variables_.ClearAndResize(num_cols); @@ -40,24 +40,56 @@ void VariablesInfo::Initialize() { void VariablesInfo::MakeBoxedVariableRelevant(bool value) { if (value == boxed_variables_are_relevant_) return; - boxed_variables_are_relevant_ = value; - for (ColIndex col(0); col < is_relevant_.size(); ++col) { - SetRelevance(col, variable_status_[col]); + if (value) { + boxed_variables_are_relevant_ = true; + for (const ColIndex col : non_basic_boxed_variables_) { + SetRelevance(col, variable_status_[col] != VariableStatus::FIXED_VALUE); + } + } else { + boxed_variables_are_relevant_ = false; + for (const ColIndex col : non_basic_boxed_variables_) { + SetRelevance(col, false); + } } } void VariablesInfo::Update(ColIndex col, VariableStatus status) { + if (status == VariableStatus::BASIC) { + UpdateToBasicStatus(col); + } else { + UpdateToNonBasicStatus(col, status); + } +} + +void VariablesInfo::UpdateToBasicStatus(ColIndex col) { variable_type_[col] = ComputeVariableType(col); + variable_status_[col] = VariableStatus::BASIC; + is_basic_.Set(col, true); + not_basic_.Set(col, false); + can_increase_.Set(col, false); + can_decrease_.Set(col, false); + non_basic_boxed_variables_.Set(col, false); + SetRelevance(col, false); +} + +void VariablesInfo::UpdateToNonBasicStatus(ColIndex col, + VariableStatus status) { + DCHECK_NE(status, VariableStatus::BASIC); + const VariableType type = ComputeVariableType(col); + variable_type_[col] = type; variable_status_[col] = status; - is_basic_.Set(col, status == VariableStatus::BASIC); - not_basic_.Set(col, status != VariableStatus::BASIC); + is_basic_.Set(col, false); + not_basic_.Set(col, true); can_increase_.Set(col, status == VariableStatus::AT_LOWER_BOUND || - status == VariableStatus::FREE); + status == VariableStatus::FREE); can_decrease_.Set(col, status == VariableStatus::AT_UPPER_BOUND || - status == VariableStatus::FREE); - SetRelevance(col, status); - non_basic_boxed_variables_.Set(col, status != VariableStatus::BASIC && - variable_type_[col] == VariableType::UPPER_AND_LOWER_BOUNDED); + status == VariableStatus::FREE); + non_basic_boxed_variables_.Set(col, + type == VariableType::UPPER_AND_LOWER_BOUNDED); + const bool relevance = status != VariableStatus::FIXED_VALUE && + (boxed_variables_are_relevant_ || + type != VariableType::UPPER_AND_LOWER_BOUNDED); + SetRelevance(col, relevance); } const VariableTypeRow& VariablesInfo::GetTypeRow() const { @@ -77,7 +109,7 @@ const DenseBitRow& VariablesInfo::GetCanDecreaseBitRow() const { } const DenseBitRow& VariablesInfo::GetIsRelevantBitRow() const { - return is_relevant_; + return relevance_; } const DenseBitRow& VariablesInfo::GetIsBasicBitRow() const { return is_basic_; } @@ -109,17 +141,13 @@ VariableType VariablesInfo::ComputeVariableType(ColIndex col) const { } } -void VariablesInfo::SetRelevance(ColIndex col, VariableStatus status) { - const bool is_relevant = (status != VariableStatus::BASIC && - status != VariableStatus::FIXED_VALUE && - (boxed_variables_are_relevant_ || - variable_type_[col] != VariableType::UPPER_AND_LOWER_BOUNDED)); - if (is_relevant_.IsSet(col) == is_relevant) return; - if (is_relevant) { - is_relevant_.Set(col); +void VariablesInfo::SetRelevance(ColIndex col, bool relevance) { + if (relevance_.IsSet(col) == relevance) return; + if (relevance) { + relevance_.Set(col); num_entries_in_relevant_columns_ += matrix_.ColumnNumEntries(col); } else { - is_relevant_.Clear(col); + relevance_.Clear(col); num_entries_in_relevant_columns_ -= matrix_.ColumnNumEntries(col); } } diff --git a/src/glop/variables_info.h b/src/glop/variables_info.h index da242342aa..432b1bcf5f 100644 --- a/src/glop/variables_info.h +++ b/src/glop/variables_info.h @@ -41,6 +41,10 @@ class VariablesInfo { // to call this if the status or the bound of a variable didn't change. void Update(ColIndex col, VariableStatus status); + // Slighlty optimized version of Update() above for the two separate cases. + void UpdateToBasicStatus(ColIndex col); + void UpdateToNonBasicStatus(ColIndex col, VariableStatus status); + // Various getter, see the corresponding member declaration below for more // information. const VariableTypeRow& GetTypeRow() const; @@ -76,7 +80,7 @@ class VariablesInfo { VariableType ComputeVariableType(ColIndex col) const; // Sets the column relevance and updates num_entries_in_relevant_columns_. - void SetRelevance(ColIndex col, VariableStatus status); + void SetRelevance(ColIndex col, bool relevance); // Problem data that should be updated from outside. const CompactSparseMatrix& matrix_; @@ -98,7 +102,7 @@ class VariablesInfo { // Indicates if we should consider this variable for entering the basis during // the simplex algorithm. Only non-fixed and non-basic columns are relevant. - DenseBitRow is_relevant_; + DenseBitRow relevance_; // Indicates if a variable is BASIC or not. There are currently two members // because the DenseBitRow class only supports a nice range-based iteration on @@ -109,7 +113,7 @@ class VariablesInfo { // Set of boxed variables that are non-basic. DenseBitRow non_basic_boxed_variables_; - // Number of entries for the relevant matrix columns (see is_relevant_). + // Number of entries for the relevant matrix columns (see relevance_). EntryIndex num_entries_in_relevant_columns_; // Whether or not a boxed variable should be considered relevant.