From c3d7946461e2dd0d3783fc4fb150c372ad67adf6 Mon Sep 17 00:00:00 2001
From: Laurent Perron <lperron@google.com>
Date: Wed, 4 Sep 2019 15:31:46 +0200
Subject: [PATCH] more work on glop

---
 ortools/glop/lu_factorization.cc |  7 ++++---
 ortools/glop/markowitz.cc        |  2 +-
 ortools/lp_data/sparse.cc        | 16 ++++++++++++----
 ortools/lp_data/sparse.h         |  4 +++-
 4 files changed, 20 insertions(+), 9 deletions(-)

diff --git a/ortools/glop/lu_factorization.cc b/ortools/glop/lu_factorization.cc
index 80bb4779ac..b4bf50de8b 100644
--- a/ortools/glop/lu_factorization.cc
+++ b/ortools/glop/lu_factorization.cc
@@ -307,12 +307,13 @@ void LuFactorization::RightSolveUWithNonZeros(ScatteredColumn* x) const {
   // If non-zeros is non-empty, we use an hypersparse solve. Note that if
   // non_zeros starts to be too big, we clear it and thus switch back to a
   // normal sparse solve.
-  upper_.ComputeRowsToConsiderInSortedOrder(&x->non_zeros);
+  upper_.ComputeRowsToConsiderInSortedOrder(&x->non_zeros, 0.1, 0.2);
   x->non_zeros_are_sorted = true;
   if (x->non_zeros.empty()) {
-    upper_.UpperSolve(&x->values);
+    transpose_upper_.TransposeLowerSolve(&x->values);
   } else {
-    upper_.HyperSparseSolveWithReversedNonZeros(&x->values, &x->non_zeros);
+    transpose_upper_.TransposeHyperSparseSolveWithReversedNonZeros(
+        &x->values, &x->non_zeros);
   }
 }
 
diff --git a/ortools/glop/markowitz.cc b/ortools/glop/markowitz.cc
index 37e66ab37a..1263fba04e 100644
--- a/ortools/glop/markowitz.cc
+++ b/ortools/glop/markowitz.cc
@@ -683,7 +683,7 @@ void MatrixNonZeroPattern::Update(RowIndex pivot_row, ColIndex pivot_col,
 
     // If the row is fully dense, there is nothing to do (the merge below will
     // not change anything). This is a small price to pay for a huge gain when
-    // the matrix become dense.
+    // the matrix becomes dense.
     if (e.coefficient() == 0.0 || row_degree_[row] == max_row_degree) continue;
     DCHECK_LT(row_degree_[row], max_row_degree);
 
diff --git a/ortools/lp_data/sparse.cc b/ortools/lp_data/sparse.cc
index 3c4ef0e32f..f7d859a6cc 100644
--- a/ortools/lp_data/sparse.cc
+++ b/ortools/lp_data/sparse.cc
@@ -1050,7 +1050,7 @@ void TriangularMatrix::PermutedLowerSparseSolve(const ColumnView& rhs,
 
   // Copy rhs into initially_all_zero_scratchpad_.
   initially_all_zero_scratchpad_.resize(num_rows_, 0.0);
-  for (SparseColumn::Entry e : rhs) {
+  for (const auto e : rhs) {
     initially_all_zero_scratchpad_[e.row()] = e.coefficient();
   }
 
@@ -1074,9 +1074,8 @@ void TriangularMatrix::PermutedLowerSparseSolve(const ColumnView& rhs,
     DCHECK_GE(row_as_col, 0);
     upper_column->SetCoefficient(permuted_row, pivot);
     DCHECK_EQ(diagonal_coefficients_[row_as_col], 1.0);
-    for (const EntryIndex i : Column(row_as_col)) {
-      initially_all_zero_scratchpad_[EntryRow(i)] -=
-          EntryCoefficient(i) * pivot;
+    for (const auto e : column(row_as_col)) {
+      initially_all_zero_scratchpad_[e.row()] -= e.coefficient() * pivot;
     }
   }
 
@@ -1306,6 +1305,15 @@ void TriangularMatrix::ComputeRowsToConsiderWithDfs(
 
 void TriangularMatrix::ComputeRowsToConsiderInSortedOrder(
     RowIndexVector* non_zero_rows) const {
+  static const Fractional kDefaultSparsityRatio = 0.025;
+  static const Fractional kDefaultNumOpsRatio = 0.05;
+  ComputeRowsToConsiderInSortedOrder(non_zero_rows, kDefaultSparsityRatio,
+                                     kDefaultNumOpsRatio);
+}
+
+void TriangularMatrix::ComputeRowsToConsiderInSortedOrder(
+    RowIndexVector* non_zero_rows, Fractional sparsity_ratio,
+    Fractional num_ops_ratio) const {
   if (non_zero_rows->empty()) return;
 
   // TODO(user): Investigate the best thresholds.
diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h
index 468f8d9511..3d54e43444 100644
--- a/ortools/lp_data/sparse.h
+++ b/ortools/lp_data/sparse.h
@@ -651,8 +651,10 @@ class TriangularMatrix : private CompactSparseMatrix {
   // sorted by rows. It is up to the client to call the direct or reverse
   // hyper-sparse solve function depending if the matrix is upper or lower
   // triangular.
+  void ComputeRowsToConsiderInSortedOrder(RowIndexVector* non_zero_rows,
+                                          Fractional sparsity_ratio,
+                                          Fractional num_ops_ratio) const;
   void ComputeRowsToConsiderInSortedOrder(RowIndexVector* non_zero_rows) const;
-
   // This is currently only used for testing. It achieves the same result as
   // PermutedLowerSparseSolve() below, but the latter exploits the sparsity of
   // rhs and is thus faster for our use case.