speedup glop; better respect the time limit

Author: lperron

ortools/glop/BUILD.bazel

@@ -286,6 +286,7 @@ cc_library(
         "//ortools/lp_data:lp_utils",
         "//ortools/lp_data:scattered_vector",
         "//ortools/util:stats",
+        "//ortools/util:time_limit",
     ],
 )
 

ortools/glop/basis_representation.h

@@ -310,6 +310,10 @@ class BasisFactorization {
   // Returns the number of updates since last refactorization.
   int NumUpdates() const { return num_updates_; }
 
+  EntryIndex NumberOfEntriesInLU() const {
+    return lu_factorization_.NumberOfEntries();
+  }
+
  private:
   // Called by ForceRefactorization() or Refactorize() or Initialize().
   Status ComputeFactorization();

ortools/glop/dual_edge_norms.cc

@@ -16,6 +16,7 @@
 #include <cstdlib>
 
 #include "ortools/lp_data/lp_utils.h"
+#include "ortools/lp_data/permutation.h"
 
 namespace operations_research {
 namespace glop {
@@ -42,8 +43,8 @@ DenseColumn::ConstView DualEdgeNorms::GetEdgeSquaredNorms() {
 void DualEdgeNorms::UpdateDataOnBasisPermutation(
     const ColumnPermutation& col_perm) {
   if (recompute_edge_squared_norms_) return;
-  ApplyColumnPermutationToRowIndexedVector(col_perm, &edge_squared_norms_,
-                                           &tmp_edge_squared_norms_);
+  ApplyColumnPermutationToRowIndexedVector(
+      col_perm.const_view(), &edge_squared_norms_, &tmp_edge_squared_norms_);
 }
 
 bool DualEdgeNorms::TestPrecision(RowIndex leaving_row,

ortools/glop/lu_factorization.cc

@@ -404,21 +404,23 @@ bool LuFactorization::LeftSolveLWithNonZeros(
   // use a different algorithm.
   ClearAndResizeVectorWithNonZeros(x->size(), result_before_permutation);
   x->swap(result_before_permutation->values);
+  const auto input = result_before_permutation->values.const_view();
+  const auto inverse_row_perm = inverse_row_perm_.const_view();
   if (nz->empty()) {
-    const RowIndex num_rows = inverse_row_perm_.size();
+    const RowIndex num_rows = inverse_row_perm.size();
     for (RowIndex row(0); row < num_rows; ++row) {
-      const Fractional value = (*result_before_permutation)[row];
+      const Fractional value = input[row];
       if (value != 0.0) {
-        const RowIndex permuted_row = inverse_row_perm_[row];
+        const RowIndex permuted_row = inverse_row_perm[row];
         (*x)[permuted_row] = value;
       }
     }
   } else {
     nz->swap(result_before_permutation->non_zeros);
     nz->reserve(result_before_permutation->non_zeros.size());
     for (const RowIndex row : result_before_permutation->non_zeros) {
-      const Fractional value = (*result_before_permutation)[row];
-      const RowIndex permuted_row = inverse_row_perm_[row];
+      const Fractional value = input[row];
+      const RowIndex permuted_row = inverse_row_perm[row];
       (*x)[permuted_row] = value;
       nz->push_back(permuted_row);
     }

ortools/glop/markowitz.cc

@@ -63,7 +63,8 @@ Status Markowitz::ComputeRowAndColumnPermutation(
   // Initialize residual_matrix_non_zero_ with the submatrix left after we
   // removed the singleton and residual singleton columns.
   residual_matrix_non_zero_.InitializeFromMatrixSubset(
-      basis_matrix, *row_perm, *col_perm, &singleton_column_, &singleton_row_);
+      basis_matrix, row_perm->const_view(), col_perm->const_view(),
+      &singleton_column_, &singleton_row_);
 
   // Perform Gaussian elimination.
   const int end_index = std::min(num_rows.value(), num_cols.value());
@@ -216,9 +217,9 @@ void Markowitz::ExtractSingletonColumns(
                                           basis_matrix.num_rows().value());
 }
 
-bool Markowitz::IsResidualSingletonColumn(const ColumnView& column,
-                                          const RowPermutation& row_perm,
-                                          RowIndex* row) {
+bool Markowitz::IsResidualSingletonColumn(
+    const ColumnView& column, StrictITISpan<RowIndex, const RowIndex> row_perm,
+    RowIndex* row) {
   int residual_degree = 0;
   for (const auto e : column) {
     if (row_perm[e.row()] != kInvalidRow) continue;
@@ -238,7 +239,9 @@ void Markowitz::ExtractResidualSingletonColumns(
   for (ColIndex col(0); col < num_cols; ++col) {
     if ((*col_perm)[col] != kInvalidCol) continue;
     const ColumnView column = basis_matrix.column(col);
-    if (!IsResidualSingletonColumn(column, *row_perm, &row)) continue;
+    if (!IsResidualSingletonColumn(column, row_perm->const_view(), &row)) {
+      continue;
+    }
     (*col_perm)[col] = ColIndex(*index);
     (*row_perm)[row] = RowIndex(*index);
     lower_.AddDiagonalOnlyColumn(1.0);
@@ -569,8 +572,10 @@ void MatrixNonZeroPattern::Reset(RowIndex num_rows, ColIndex num_cols) {
 }
 
 void MatrixNonZeroPattern::InitializeFromMatrixSubset(
-    const CompactSparseMatrixView& basis_matrix, const RowPermutation& row_perm,
-    const ColumnPermutation& col_perm, std::vector<ColIndex>* singleton_columns,
+    const CompactSparseMatrixView& basis_matrix,
+    StrictITISpan<RowIndex, const RowIndex> row_perm,
+    StrictITISpan<ColIndex, const ColIndex> col_perm,
+    std::vector<ColIndex>* singleton_columns,
     std::vector<RowIndex>* singleton_rows) {
   const ColIndex num_cols = basis_matrix.num_cols();
   const RowIndex num_rows = basis_matrix.num_rows();

ortools/glop/markowitz.h

@@ -116,11 +116,12 @@ class MatrixNonZeroPattern {
   // Resets the pattern to the one of the given matrix but only for the
   // rows/columns whose given permutation is kInvalidRow or kInvalidCol.
   // This also fills the singleton columns/rows with the corresponding entries.
-  void InitializeFromMatrixSubset(const CompactSparseMatrixView& basis_matrix,
-                                  const RowPermutation& row_perm,
-                                  const ColumnPermutation& col_perm,
-                                  std::vector<ColIndex>* singleton_columns,
-                                  std::vector<RowIndex>* singleton_rows);
+  void InitializeFromMatrixSubset(
+      const CompactSparseMatrixView& basis_matrix,
+      StrictITISpan<RowIndex, const RowIndex> row_perm,
+      StrictITISpan<ColIndex, const ColIndex> col_perm,
+      std::vector<ColIndex>* singleton_columns,
+      std::vector<RowIndex>* singleton_rows);
 
   // Adds a non-zero entry to the matrix. There should be no duplicates.
   void AddEntry(RowIndex row, ColIndex col);
@@ -377,8 +378,9 @@ class Markowitz {
 
   // 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);
+  bool IsResidualSingletonColumn(
+      const ColumnView& column,
+      StrictITISpan<RowIndex, const RowIndex> 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

ortools/glop/primal_edge_norms.cc

@@ -156,16 +156,27 @@ void PrimalEdgeNorms::ComputeMatrixColumnNorms() {
 void PrimalEdgeNorms::ComputeEdgeSquaredNorms() {
   SCOPED_TIME_STAT(&stats_);
 
+  // time_limit_->LimitReached() can be costly sometimes, so we only do that
+  // if we feel this will be slow anyway.
+  const bool test_limit = (time_limit_ != nullptr) &&
+                          basis_factorization_.NumberOfEntriesInLU() > 10'000;
+
   // Since we will do a lot of inversions, it is better to be as efficient and
   // precise as possible by refactorizing the basis.
   DCHECK(basis_factorization_.IsRefactorized());
-  edge_squared_norms_.resize(compact_matrix_.num_cols(), 0.0);
+  edge_squared_norms_.resize(compact_matrix_.num_cols(), 1.0);
   for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
     // Note the +1.0 in the squared norm for the component of the edge on the
     // 'entering_col'.
     edge_squared_norms_[col] = 1.0 + basis_factorization_.RightSolveSquaredNorm(
                                          compact_matrix_.column(col));
+
+    // This operation can be costly, and we abort if we are stuck here.
+    // Note that we still mark edges as "recomputed" otherwise we can runs into
+    // some DCHECK before we actually abort the solve.
+    if (test_limit && time_limit_->LimitReached()) break;
   }
+
   recompute_edge_squared_norms_ = false;
 }
 

ortools/glop/primal_edge_norms.h

@@ -27,6 +27,7 @@
 #include "ortools/lp_data/scattered_vector.h"
 #include "ortools/lp_data/sparse.h"
 #include "ortools/util/stats.h"
+#include "ortools/util/time_limit.h"
 
 namespace operations_research {
 namespace glop {
@@ -146,6 +147,8 @@ class PrimalEdgeNorms {
     return DeterministicTimeForFpOperations(num_operations_);
   }
 
+  void SetTimeLimit(TimeLimit* time_limit) { time_limit_ = time_limit; }
+
  private:
   // Statistics about this class.
   struct Stats : public StatsGroup {
@@ -192,6 +195,7 @@ class PrimalEdgeNorms {
   const CompactSparseMatrix& compact_matrix_;
   const VariablesInfo& variables_info_;
   const BasisFactorization& basis_factorization_;
+  TimeLimit* time_limit_ = nullptr;
 
   // Internal data.
   GlopParameters parameters_;

ortools/glop/rank_one_update.h

@@ -177,7 +177,7 @@ class RankOneUpdateFactorization {
     }
 
     // y->is_non_zero is always all false before and after this code.
-    DCHECK(IsAllFalse(y->is_non_zero));
+    DCHECK(y->is_non_zero.IsAllFalse());
     y->RepopulateSparseMask();
     bool use_dense = y->ShouldUseDenseIteration(hypersparse_ratio_);
     for (int i = elementary_matrices_.size() - 1; i >= 0; --i) {
@@ -213,7 +213,7 @@ class RankOneUpdateFactorization {
     }
 
     // d->is_non_zero is always all false before and after this code.
-    DCHECK(IsAllFalse(d->is_non_zero));
+    DCHECK(d->is_non_zero.IsAllFalse());
     d->RepopulateSparseMask();
     bool use_dense = d->ShouldUseDenseIteration(hypersparse_ratio_);
     const size_t end = elementary_matrices_.size();

ortools/glop/revised_simplex.cc

@@ -220,6 +220,7 @@ ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal(
       [this, time_limit]() { AdvanceDeterministicTime(time_limit); });
 
   SOLVER_LOG(logger_, "");
+  primal_edge_norms_.SetTimeLimit(time_limit);
   if (logger_->LoggingIsEnabled()) {
     DisplayBasicVariableStatistics();
   }
@@ -2563,13 +2564,15 @@ void RevisedSimplex::PermuteBasis() {
   if (col_perm.empty()) return;
 
   // Permute basis_.
-  ApplyColumnPermutationToRowIndexedVector(col_perm, &basis_, &tmp_basis_);
+  ApplyColumnPermutationToRowIndexedVector(col_perm.const_view(), &basis_,
+                                           &tmp_basis_);
 
   // Permute dual_pricing_vector_ if needed.
   if (!dual_pricing_vector_.empty()) {
     // TODO(user): We need to permute dual_prices_ too now, we recompute
     // everything one each basis factorization, so this don't matter.
-    ApplyColumnPermutationToRowIndexedVector(col_perm, &dual_pricing_vector_,
+    ApplyColumnPermutationToRowIndexedVector(col_perm.const_view(),
+                                             &dual_pricing_vector_,
                                              &tmp_dual_pricing_vector_);
   }
 

ortools/glop/update_row.cc

@@ -106,20 +106,21 @@ void UpdateRow::ComputeUpdateRow(RowIndex leaving_row) {
     // small entries (no complexity changes).
     const Fractional drop_tolerance = parameters_.drop_tolerance();
     unit_row_left_inverse_filtered_non_zeros_.clear();
-    const auto view = transposed_matrix_.view();
+    const auto matrix_view = transposed_matrix_.view();
     if (unit_row_left_inverse_.non_zeros.empty()) {
       const ColIndex size = unit_row_left_inverse_.values.size();
+      const auto values = unit_row_left_inverse_.values.view();
       for (ColIndex col(0); col < size; ++col) {
-        if (std::abs(unit_row_left_inverse_.values[col]) > drop_tolerance) {
+        if (std::abs(values[col]) > drop_tolerance) {
           unit_row_left_inverse_filtered_non_zeros_.push_back(col);
-          num_row_wise_entries += view.ColumnNumEntries(col);
+          num_row_wise_entries += matrix_view.ColumnNumEntries(col);
         }
       }
     } else {
       for (const auto e : unit_row_left_inverse_) {
         if (std::abs(e.coefficient()) > drop_tolerance) {
           unit_row_left_inverse_filtered_non_zeros_.push_back(e.column());
-          num_row_wise_entries += view.ColumnNumEntries(e.column());
+          num_row_wise_entries += matrix_view.ColumnNumEntries(e.column());
         }
       }
     }

ortools/lp_data/lp_utils.h

@@ -249,12 +249,17 @@ inline void PermuteWithScratchpad(
   const IndexType size = input_output->size();
   zero_scratchpad->swap(*input_output);
   input_output->resize(size, 0.0);
-  for (IndexType index(0); index < size; ++index) {
-    const Fractional value = (*zero_scratchpad)[index];
+
+  // Caching the pointers help.
+  const Fractional* const input = zero_scratchpad->data();
+  const IndexType* const perm = permutation.data();
+  Fractional* const output = input_output->data();
+  for (int i = 0; i < size; ++i) {
+    DCHECK_GE(perm[i], 0);
+    DCHECK_LT(perm[i], size);
+    const Fractional value = input[i];
     if (value != 0.0) {
-      const IndexType permuted_index(
-          permutation[PermutationIndexType(index.value())].value());
-      (*input_output)[permuted_index] = value;
+      output[perm[i].value()] = value;
     }
   }
   zero_scratchpad->assign(size, 0.0);