From db2d441ca2f711b86a7bad468f710eb9219888f4 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 20 Aug 2019 15:33:47 -0700 Subject: [PATCH] more work on GLOP internals --- ortools/glop/lu_factorization.cc | 8 ++--- ortools/glop/markowitz.cc | 51 ++++++++++++++++++-------------- ortools/glop/markowitz.h | 5 ++++ ortools/glop/preprocessor.cc | 3 +- ortools/glop/preprocessor.h | 1 + ortools/glop/revised_simplex.h | 4 +-- ortools/glop/variable_values.cc | 2 +- ortools/glop/variables_info.h | 2 +- ortools/lp_data/sparse.cc | 37 +++++++++++++---------- ortools/lp_data/sparse.h | 15 +++++++--- 10 files changed, 77 insertions(+), 51 deletions(-) diff --git a/ortools/glop/lu_factorization.cc b/ortools/glop/lu_factorization.cc index dd479529df..341206bff4 100644 --- a/ortools/glop/lu_factorization.cc +++ b/ortools/glop/lu_factorization.cc @@ -30,10 +30,10 @@ LuFactorization::LuFactorization() void LuFactorization::Clear() { SCOPED_TIME_STAT(&stats_); - lower_.Reset(RowIndex(0)); - upper_.Reset(RowIndex(0)); - transpose_upper_.Reset(RowIndex(0)); - transpose_lower_.Reset(RowIndex(0)); + lower_.Reset(RowIndex(0), ColIndex(0)); + upper_.Reset(RowIndex(0), ColIndex(0)); + transpose_upper_.Reset(RowIndex(0), ColIndex(0)); + transpose_lower_.Reset(RowIndex(0), ColIndex(0)); is_identity_factorization_ = true; col_perm_.clear(); row_perm_.clear(); diff --git a/ortools/glop/markowitz.cc b/ortools/glop/markowitz.cc index 7a6d4b0173..37e66ab37a 100644 --- a/ortools/glop/markowitz.cc +++ b/ortools/glop/markowitz.cc @@ -16,6 +16,7 @@ #include #include "absl/strings/str_format.h" +#include "ortools/lp_data/lp_types.h" #include "ortools/lp_data/lp_utils.h" #include "ortools/lp_data/sparse.h" @@ -37,8 +38,8 @@ Status Markowitz::ComputeRowAndColumnPermutation( basis_matrix_ = &basis_matrix; // Initialize all the matrices. - lower_.Reset(num_rows); - upper_.Reset(num_rows); + lower_.Reset(num_rows, num_cols); + upper_.Reset(num_rows, num_cols); permuted_lower_.Reset(num_cols); permuted_upper_.Reset(num_cols); permuted_lower_column_needs_solve_.assign(num_cols, false); @@ -213,30 +214,34 @@ void Markowitz::ExtractSingletonColumns( num_cols.value()); } +bool Markowitz::IsResidualSingletonColumn(const ColumnView& column, + const RowPermutation& row_perm, + RowIndex* row) { + int residual_degree = 0; + for (const auto e : column) { + if (row_perm[e.row()] != kInvalidRow) continue; + ++residual_degree; + if (residual_degree > 1) return false; + *row = e.row(); + } + return residual_degree == 1; +} + void Markowitz::ExtractResidualSingletonColumns( const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm, ColumnPermutation* col_perm, int* index) { SCOPED_TIME_STAT(&stats_); const ColIndex num_cols = basis_matrix.num_cols(); + RowIndex row = kInvalidRow; for (ColIndex col(0); col < num_cols; ++col) { if ((*col_perm)[col] != kInvalidCol) continue; const ColumnView& column = basis_matrix.column(col); - int residual_degree = 0; - RowIndex row; - for (const SparseColumn::Entry e : column) { - if ((*row_perm)[e.row()] == kInvalidRow) { - ++residual_degree; - if (residual_degree > 1) break; - row = e.row(); - } - } - if (residual_degree == 1) { - (*col_perm)[col] = ColIndex(*index); - (*row_perm)[row] = RowIndex(*index); - lower_.AddDiagonalOnlyColumn(1.0); - upper_.AddTriangularColumn(column, row); - ++(*index); - } + if (!IsResidualSingletonColumn(column, *row_perm, &row)) continue; + (*col_perm)[col] = ColIndex(*index); + (*row_perm)[row] = RowIndex(*index); + lower_.AddDiagonalOnlyColumn(1.0); + upper_.AddTriangularColumn(column, row); + ++(*index); } stats_.basis_residual_singleton_column_ratio.Add(static_cast(*index) / num_cols.value()); @@ -543,12 +548,12 @@ void MatrixNonZeroPattern::Clear() { } void MatrixNonZeroPattern::Reset(RowIndex num_rows, ColIndex num_cols) { - Clear(); - row_degree_.resize(num_rows, 0); - col_degree_.resize(num_cols, 0); + row_degree_.AssignToZero(num_rows); + col_degree_.AssignToZero(num_cols); + row_non_zero_.clear(); row_non_zero_.resize(num_rows.value()); - deleted_columns_.resize(num_cols, false); - bool_scratchpad_.resize(num_cols, false); + deleted_columns_.assign(num_cols, false); + bool_scratchpad_.assign(num_cols, false); num_non_deleted_columns_ = num_cols; } diff --git a/ortools/glop/markowitz.h b/ortools/glop/markowitz.h index 8143133a5b..74c72f0205 100644 --- a/ortools/glop/markowitz.h +++ b/ortools/glop/markowitz.h @@ -347,6 +347,11 @@ class Markowitz { const CompactSparseMatrixView& basis_matrix, RowPermutation* row_perm, ColumnPermutation* col_perm, int* index); + // Helper function for determining if a column is a residual singleton column. + // If it is, RowIndex* row contains the index of the single residual edge. + bool IsResidualSingletonColumn(const ColumnView& column, + const RowPermutation& row_perm, RowIndex* row); + // Returns the column of the current residual matrix with an index 'col' in // the initial matrix. We compute it by solving a linear system with the // current lower_ and the last computed column 'col' of a previous residual diff --git a/ortools/glop/preprocessor.cc b/ortools/glop/preprocessor.cc index 240840bcac..6a9d7a71e5 100644 --- a/ortools/glop/preprocessor.cc +++ b/ortools/glop/preprocessor.cc @@ -44,7 +44,8 @@ Preprocessor::Preprocessor(const GlopParameters* parameters) : status_(ProblemStatus::INIT), parameters_(*parameters), in_mip_context_(false), - time_limit_(TimeLimit::Infinite().get()) {} + infinite_time_limit_(TimeLimit::Infinite()), + time_limit_(infinite_time_limit_.get()) {} Preprocessor::~Preprocessor() {} // -------------------------------------------------------- diff --git a/ortools/glop/preprocessor.h b/ortools/glop/preprocessor.h index 7323279760..3a5b307b16 100644 --- a/ortools/glop/preprocessor.h +++ b/ortools/glop/preprocessor.h @@ -90,6 +90,7 @@ class Preprocessor { ProblemStatus status_; const GlopParameters& parameters_; bool in_mip_context_; + std::unique_ptr infinite_time_limit_; TimeLimit* time_limit_; }; diff --git a/ortools/glop/revised_simplex.h b/ortools/glop/revised_simplex.h index b36bd20cca..3afc2cb664 100644 --- a/ortools/glop/revised_simplex.h +++ b/ortools/glop/revised_simplex.h @@ -247,8 +247,8 @@ class RevisedSimplex { // Propagates parameters_ to all the other classes that need it. // // TODO(user): Maybe a better design is for them to have a reference to a - // unique parameters object? it will clutter a bit more these classes - // contructor though. + // unique parameters object? It will clutter a bit more these classes' + // constructor though. void PropagateParameters(); // Returns a std::string containing the same information as with diff --git a/ortools/glop/variable_values.cc b/ortools/glop/variable_values.cc index 5651803303..dc8a6090e3 100644 --- a/ortools/glop/variable_values.cc +++ b/ortools/glop/variable_values.cc @@ -150,7 +150,7 @@ void VariableValues::UpdateOnPivoting(const ScatteredColumn& direction, // - The rows for which direction_[row] < tolerance. // - The non-zeros of direction_ignored_position_ in case of degeneracy. // Such positions may result in basic variables going out of their bounds by - // more than the allowed tolerance. We could choose not to update theses + // more than the allowed tolerance. We could choose not to update these // variables or not make them take out-of-bound values, but this would // introduce artificial errors. diff --git a/ortools/glop/variables_info.h b/ortools/glop/variables_info.h index bfc4c31389..4cc138acea 100644 --- a/ortools/glop/variables_info.h +++ b/ortools/glop/variables_info.h @@ -43,7 +43,7 @@ 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. + // Slightly optimized version of Update() above for the two separate cases. void UpdateToBasicStatus(ColIndex col); void UpdateToNonBasicStatus(ColIndex col, VariableStatus status); diff --git a/ortools/lp_data/sparse.cc b/ortools/lp_data/sparse.cc index dfcdc9197f..3c4ef0e32f 100644 --- a/ortools/lp_data/sparse.cc +++ b/ortools/lp_data/sparse.cc @@ -16,6 +16,7 @@ #include #include "absl/strings/str_format.h" +#include "ortools/base/logging.h" #include "ortools/lp_data/lp_types.h" #include "ortools/lp_data/permutation.h" @@ -520,12 +521,16 @@ void CompactSparseMatrix::Reset(RowIndex num_rows) { starts_.push_back(EntryIndex(0)); } -void TriangularMatrix::Reset(RowIndex num_rows) { +void TriangularMatrix::Reset(RowIndex num_rows, ColIndex col_capacity) { CompactSparseMatrix::Reset(num_rows); first_non_identity_column_ = 0; - diagonal_coefficients_.clear(); all_diagonal_coefficients_are_one_ = true; - pruned_ends_.clear(); + + pruned_ends_.resize(col_capacity); + diagonal_coefficients_.resize(col_capacity); + starts_.resize(col_capacity + 1); + // Non-zero entries in the first column always have an offset of 0. + starts_[ColIndex(0)] = 0; } ColIndex CompactSparseMatrix::AddDenseColumn(const DenseColumn& dense_column) { @@ -608,22 +613,24 @@ Fractional CompactSparseMatrixView::ComputeInfinityNorm() const { // matrix for the next column addition. void TriangularMatrix::CloseCurrentColumn(Fractional diagonal_value) { DCHECK_NE(diagonal_value, 0.0); + // The vectors diagonal_coefficients, pruned_ends, and starts_ should have all + // been preallocated by a call to SetTotalNumberOfColumns(). + DCHECK_LT(num_cols_, diagonal_coefficients_.size()); + diagonal_coefficients_[num_cols_] = diagonal_value; + + // TODO(user): This is currently not used by all matrices. It will be good + // to fill it only when needed. + DCHECK_LT(num_cols_, pruned_ends_.size()); + pruned_ends_[num_cols_] = coefficients_.size(); ++num_cols_; - starts_.push_back(coefficients_.size()); - diagonal_coefficients_.push_back(diagonal_value); - DCHECK_EQ(num_cols_, diagonal_coefficients_.size()); - DCHECK_EQ(num_cols_ + 1, starts_.size()); + DCHECK_LT(num_cols_, starts_.size()); + starts_[num_cols_] = coefficients_.size(); if (first_non_identity_column_ == num_cols_ - 1 && coefficients_.empty() && diagonal_value == 1.0) { first_non_identity_column_ = num_cols_; } - if (all_diagonal_coefficients_are_one_) { - all_diagonal_coefficients_are_one_ = (diagonal_value == 1.0); - } - - // TODO(user): This is currently not used by all matrices. It will be good - // to fill it only when needed. - pruned_ends_.push_back(coefficients_.size()); + all_diagonal_coefficients_are_one_ = + all_diagonal_coefficients_are_one_ && (diagonal_value == 1.0); } void TriangularMatrix::AddDiagonalOnlyColumn(Fractional diagonal_value) { @@ -675,7 +682,7 @@ void TriangularMatrix::AddTriangularColumnWithGivenDiagonalEntry( void TriangularMatrix::PopulateFromTriangularSparseMatrix( const SparseMatrix& input) { - Reset(input.num_rows()); + Reset(input.num_rows(), input.num_cols()); for (ColIndex col(0); col < input.num_cols(); ++col) { AddTriangularColumn(ColumnView(input.column(col)), ColToRowIndex(col)); } diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h index bbff3d33a7..468f8d9511 100644 --- a/ortools/lp_data/sparse.h +++ b/ortools/lp_data/sparse.h @@ -506,7 +506,6 @@ class TriangularMatrix : private CompactSparseMatrix { // Only a subset of the functions from CompactSparseMatrix are exposed (note // the private inheritance). They are extended to deal with diagonal // coefficients properly. - void Reset(RowIndex num_rows); void PopulateFromTranspose(const TriangularMatrix& input); void Swap(TriangularMatrix* other); bool IsEmpty() const { return diagonal_coefficients_.empty(); } @@ -516,17 +515,25 @@ class TriangularMatrix : private CompactSparseMatrix { return EntryIndex(num_cols_.value()) + coefficients_.size(); } + // On top of the CompactSparseMatrix functionality, TriangularMatrix::Reset() + // also pre-allocates space of size col_size for a number of internal vectors. + // This helps reduce costly push_back operations for large problems. + // + // WARNING: Reset() must be called with a sufficiently large col_capacity + // prior to any Add* calls (e.g., AddTriangularColumn). + void Reset(RowIndex num_rows, ColIndex col_capacity); + // Constructs a triangular matrix from the given SparseMatrix. The input is // assumed to be lower or upper triangular without any permutations. This is // checked in debug mode. void PopulateFromTriangularSparseMatrix(const SparseMatrix& input); - // Functions to create a triangular matrix incrementally, columns by columns. - // A client need to call Reset(num_rows) first, and then each column must be + // Functions to create a triangular matrix incrementally, column by column. + // A client needs to call Reset(num_rows) first, and then each column must be // added by calling one of the 3 functions below. // // Note that the row indices of the columns are allowed to be permuted: the - // diagonal entry of the column #col not beeing necessarily on the row #col. + // diagonal entry of the column #col not being necessarily on the row #col. // This is why these functions require the 'diagonal_row' parameter. The // permutation can be fixed at the end by a call to // ApplyRowPermutationToNonDiagonalEntries() or accounted directly in the case