fix glop on 1 mittleman instance

src/glop/basis_representation.cc

@@ -185,7 +185,9 @@ BasisFactorization::BasisFactorization(const MatrixView& matrix,
       num_updates_(0),
       eta_factorization_(),
       lu_factorization_(),
-      deterministic_time_(0.0) {}
+      deterministic_time_(0.0) {
+  SetParameters(parameters_);
+}
 
 BasisFactorization::~BasisFactorization() {}
 

src/glop/initial_basis.cc

@@ -96,18 +96,18 @@ void InitialBasis::CompleteBixbyBasis(ColIndex num_cols,
   }
 }
 
-void InitialBasis::CompleteTriangularPrimalBasis(ColIndex num_cols,
+bool InitialBasis::CompleteTriangularPrimalBasis(ColIndex num_cols,
                                                  RowToColMapping* basis) {
-  CompleteTriangularBasis<false>(num_cols, basis);
+  return CompleteTriangularBasis<false>(num_cols, basis);
 }
 
-void InitialBasis::CompleteTriangularDualBasis(ColIndex num_cols,
+bool InitialBasis::CompleteTriangularDualBasis(ColIndex num_cols,
                                                RowToColMapping* basis) {
-  CompleteTriangularBasis<true>(num_cols, basis);
+  return CompleteTriangularBasis<true>(num_cols, basis);
 }
 
 template <bool only_allow_zero_cost_column>
-void InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
+bool InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
                                            RowToColMapping* basis) {
   // Initialize can_be_replaced.
   const RowIndex num_rows = matrix_.num_rows();
@@ -150,6 +150,14 @@ void InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
       queue(residual_singleton_column.begin(), residual_singleton_column.end(),
             triangular_column_comparator_);
 
+  // If the the product magnitude of the diagonal coefficients become smaller
+  // than a given threshold, we will assume that this method returns an instable
+  // first basis. The threshold is somewhat arbitrary and is mainly here to
+  // avoid an infinite inverse product which will trigger floating point
+  // exceptions in other part of the code.
+  const double kMinimumProductMagnitude = 1e-100;
+  double partial_diagonal_product = 1.0;
+
   // Process the residual singleton columns by priority and add them to the
   // basis if their "diagonal" coefficient is not too small.
   while (!queue.empty()) {
@@ -174,6 +182,14 @@ void InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
     if (fabs(coeff) < kStabilityThreshold * max_magnitude) continue;
     DCHECK_NE(kInvalidRow, row);
 
+    partial_diagonal_product *= coeff;
+    if (fabs(partial_diagonal_product) < kMinimumProductMagnitude) {
+      LOG(WARNING) << "Numerical difficulties detected. Product of "
+                   << "diagonal coefficients is currently equals to "
+                   << partial_diagonal_product;
+      break;
+    }
+
     // Use this candidate column in the basis.
     (*basis)[row] = candidate;
     can_be_replaced[row] = false;
@@ -186,6 +202,8 @@ void InitialBasis::CompleteTriangularBasis(ColIndex num_cols,
       }
     }
   }
+
+  return fabs(partial_diagonal_product) >= kMinimumProductMagnitude;
 }
 
 void InitialBasis::ComputeCandidates(ColIndex num_cols,

src/glop/initial_basis.h

@@ -59,8 +59,11 @@ class InitialBasis {
   // one. This function usually produces better initial bases. The dual version
   // just restricts the possible entering columns to the ones with a cost of 0.0
   // in order to always start with the all-zeros vector of dual values.
-  void CompleteTriangularPrimalBasis(ColIndex num_cols, RowToColMapping* basis);
-  void CompleteTriangularDualBasis(ColIndex num_cols, RowToColMapping* basis);
+  //
+  // Returns false if an error occurred during the algorithm (numerically
+  // instable basis).
+  bool CompleteTriangularPrimalBasis(ColIndex num_cols, RowToColMapping* basis);
+  bool CompleteTriangularDualBasis(ColIndex num_cols, RowToColMapping* basis);
 
   // Visible for testing. Computes a list of candidate column indices out of the
   // fist num_candidate_columns of A and sorts them using the
@@ -71,7 +74,7 @@ class InitialBasis {
  private:
   // Internal implementation of the Primal/Dual CompleteTriangularBasis().
   template <bool only_allow_zero_cost_column>
-  void CompleteTriangularBasis(ColIndex num_cols, RowToColMapping* basis);
+  bool CompleteTriangularBasis(ColIndex num_cols, RowToColMapping* basis);
 
   // Returns an integer representing the order (the lower the better)
   // between column categories (known as C2, C3 or C4 in the paper).

src/glop/revised_simplex.cc

@@ -832,7 +832,14 @@ Status RevisedSimplex::CreateInitialBasis() {
       } else if (parameters_.initial_basis() == GlopParameters::TRIANGULAR) {
         // Note the use of num_cols_ here because this algorithm
         // benefits from treating fixed slack columns like any other column.
-        initial_basis.CompleteTriangularPrimalBasis(num_cols_, &basis);
+        RowToColMapping basis_copy = basis;
+        if (!initial_basis.CompleteTriangularPrimalBasis(num_cols_, &basis)) {
+          LOG(WARNING) << "Reverting to Bixby's initial basis algorithm.";
+          basis = basis_copy;
+          if (parameters_.use_scaling()) {
+            initial_basis.CompleteBixbyBasis(first_slack_col_, &basis);
+          }
+        }
       } else {
         LOG(WARNING) << "Unsupported initial_basis parameters: "
                      << parameters_.initial_basis();