Use Fractional everywhere

Author: lperron

ortools/glop/lp_solver.cc

@@ -325,7 +325,7 @@ Fractional ProblemObjectiveValue(const LinearProgram& lp, Fractional value) {
 // Returns the allowed error magnitude for something that should evaluate to
 // value under the given tolerance.
 Fractional AllowedError(Fractional tolerance, Fractional value) {
-  return tolerance * std::max(1.0, std::abs(value));
+  return tolerance * std::max(Fractional(1.0), std::abs(value));
 }
 }  // namespace
 
@@ -492,15 +492,15 @@ bool LPSolver::IsOptimalSolutionOnFacet(const LinearProgram& lp) {
   // primal values.
   // TODO(user): investigate whether to use the tolerances defined in
   // parameters.proto.
-  const double kReducedCostTolerance = 1e-9;
-  const double kBoundTolerance = 1e-7;
+  const Fractional kReducedCostTolerance = 1e-9;
+  const Fractional kBoundTolerance = 1e-7;
   const ColIndex num_cols = lp.num_variables();
   for (ColIndex col(0); col < num_cols; ++col) {
     if (variable_statuses_[col] == VariableStatus::FIXED_VALUE) continue;
     const Fractional lower_bound = lp.variable_lower_bounds()[col];
     const Fractional upper_bound = lp.variable_upper_bounds()[col];
     const Fractional value = primal_values_[col];
-    if (AreWithinAbsoluteTolerance(reduced_costs_[col], 0.0,
+    if (AreWithinAbsoluteTolerance(reduced_costs_[col], Fractional(0.0),
                                    kReducedCostTolerance) &&
         (AreWithinAbsoluteTolerance(value, lower_bound, kBoundTolerance) ||
          AreWithinAbsoluteTolerance(value, upper_bound, kBoundTolerance))) {
@@ -513,7 +513,7 @@ bool LPSolver::IsOptimalSolutionOnFacet(const LinearProgram& lp) {
     const Fractional lower_bound = lp.constraint_lower_bounds()[row];
     const Fractional upper_bound = lp.constraint_upper_bounds()[row];
     const Fractional activity = constraint_activities_[row];
-    if (AreWithinAbsoluteTolerance(dual_values_[row], 0.0,
+    if (AreWithinAbsoluteTolerance(dual_values_[row], Fractional(0.0),
                                    kReducedCostTolerance) &&
         (AreWithinAbsoluteTolerance(activity, lower_bound, kBoundTolerance) ||
          AreWithinAbsoluteTolerance(activity, upper_bound, kBoundTolerance))) {
@@ -739,7 +739,8 @@ bool LPSolver::IsProblemSolutionConsistent(
       case VariableStatus::AT_UPPER_BOUND:
         // TODO(user): revert to an exact comparison once the bug causing this
         // to fail has been fixed.
-        if (!AreWithinAbsoluteTolerance(value, ub, 1e-7) || lb == ub) {
+        if (!AreWithinAbsoluteTolerance(value, ub, Fractional(1e-7)) ||
+            lb == ub) {
           LogVariableStatusError(col, value, status, lb, ub);
           return false;
         }
@@ -1002,7 +1003,7 @@ double LPSolver::ComputeMaxExpectedObjectiveError(const LinearProgram& lp) {
 
 double LPSolver::ComputePrimalValueInfeasibility(const LinearProgram& lp,
                                                  bool* is_too_large) {
-  double infeasibility = 0.0;
+  Fractional infeasibility = 0.0;
   const Fractional tolerance = parameters_.solution_feasibility_tolerance();
   const ColIndex num_cols = lp.num_variables();
   for (ColIndex col(0); col < num_cols; ++col) {
@@ -1032,7 +1033,7 @@ double LPSolver::ComputePrimalValueInfeasibility(const LinearProgram& lp,
 
 double LPSolver::ComputeActivityInfeasibility(const LinearProgram& lp,
                                               bool* is_too_large) {
-  double infeasibility = 0.0;
+  Fractional infeasibility = 0.0;
   int num_problematic_rows(0);
   const RowIndex num_rows = lp.num_constraints();
   const Fractional tolerance = parameters_.solution_feasibility_tolerance();
@@ -1085,7 +1086,7 @@ double LPSolver::ComputeDualValueInfeasibility(const LinearProgram& lp,
                                                bool* is_too_large) {
   const Fractional allowed_error = parameters_.solution_feasibility_tolerance();
   const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0;
-  double infeasibility = 0.0;
+  Fractional infeasibility = 0.0;
   const RowIndex num_rows = lp.num_constraints();
   for (RowIndex row(0); row < num_rows; ++row) {
     const Fractional dual_value = dual_values_[row];
@@ -1108,7 +1109,7 @@ double LPSolver::ComputeDualValueInfeasibility(const LinearProgram& lp,
 double LPSolver::ComputeReducedCostInfeasibility(const LinearProgram& lp,
                                                  bool* is_too_large) {
   const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0;
-  double infeasibility = 0.0;
+  Fractional infeasibility = 0.0;
   const ColIndex num_cols = lp.num_variables();
   const Fractional tolerance = parameters_.solution_feasibility_tolerance();
   for (ColIndex col(0); col < num_cols; ++col) {

ortools/glop/preprocessor.cc

@@ -2430,7 +2430,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable(
     const Fractional coefficient_ratio = coefficient / matrix_entry.coeff;
     // Check if coefficient_ratio is integer.
     if (!IsIntegerWithinTolerance(
-            coefficient_ratio, parameters_.solution_feasibility_tolerance())) {
+            coefficient_ratio,
+            Fractional(parameters_.solution_feasibility_tolerance()))) {
       return false;
     }
   }
@@ -2439,7 +2440,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable(
   if (IsFinite(constraint_lb)) {
     const Fractional lower_bound_ratio = constraint_lb / matrix_entry.coeff;
     if (!IsIntegerWithinTolerance(
-            lower_bound_ratio, parameters_.solution_feasibility_tolerance())) {
+            lower_bound_ratio,
+            Fractional(parameters_.solution_feasibility_tolerance()))) {
       return false;
     }
   }
@@ -2448,7 +2450,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable(
   if (IsFinite(constraint_ub)) {
     const Fractional upper_bound_ratio = constraint_ub / matrix_entry.coeff;
     if (!IsIntegerWithinTolerance(
-            upper_bound_ratio, parameters_.solution_feasibility_tolerance())) {
+            upper_bound_ratio,
+            Fractional(parameters_.solution_feasibility_tolerance()))) {
       return false;
     }
   }

ortools/glop/preprocessor.h

@@ -83,13 +83,13 @@ class Preprocessor {
   // tolerance).
   bool IsSmallerWithinFeasibilityTolerance(Fractional a, Fractional b) const {
     return ::operations_research::IsSmallerWithinTolerance(
-        a, b, parameters_.solution_feasibility_tolerance());
+        a, b, Fractional(parameters_.solution_feasibility_tolerance()));
   }
   bool IsSmallerWithinPreprocessorZeroTolerance(Fractional a,
                                                 Fractional b) const {
     // TODO(user): use an absolute tolerance here to be even more defensive?
     return ::operations_research::IsSmallerWithinTolerance(
-        a, b, parameters_.preprocessor_zero_tolerance());
+        a, b, Fractional(parameters_.preprocessor_zero_tolerance()));
   }
 
   ProblemStatus status_;

ortools/glop/primal_edge_norms.cc

@@ -238,11 +238,14 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col,
   const Fractional pivot = -direction[leaving_row];
   DCHECK_NE(pivot, 0.0);
 
+  const ColIndex first_slack =
+      compact_matrix_.num_cols() - RowToColIndex(compact_matrix_.num_rows());
+
   // Note that this should be precise because of the call to
   // TestEnteringEdgeNormPrecision().
   const Fractional entering_squared_norm = edge_squared_norms_[entering_col];
   const Fractional leaving_squared_norm =
-      std::max(1.0, entering_squared_norm / Square(pivot));
+      std::max(Fractional(1.0), entering_squared_norm / Square(pivot));
 
   int stat_lower_bounded_norms = 0;
   const Fractional factor = 2.0 / pivot;
@@ -253,7 +256,9 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col,
   for (const ColIndex col : update_row.GetNonZeroPositions()) {
     const Fractional coeff = update_row.GetCoefficient(col);
     const Fractional scalar_product =
-        view.ColumnScalarProduct(col, direction_left_inverse);
+        col >= first_slack
+            ? direction_left_inverse_[col - first_slack]
+            : view.ColumnScalarProduct(col, direction_left_inverse);
     num_operations_ += view.ColumnNumEntries(col).value();
 
     // Update the edge squared norm of this column. Note that the update
@@ -288,7 +293,7 @@ void PrimalEdgeNorms::UpdateDevexWeights(
   const Fractional entering_norm = sqrt(PreciseSquaredNorm(direction));
   const Fractional pivot_magnitude = std::abs(direction[leaving_row]);
   const Fractional leaving_norm =
-      std::max(1.0, entering_norm / pivot_magnitude);
+      std::max(Fractional(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 = std::abs(coeff) * leaving_norm;

ortools/glop/reduced_costs.cc

@@ -26,17 +26,12 @@
 #include "ortools/glop/update_row.h"
 #include "ortools/glop/variables_info.h"
 #include "ortools/lp_data/lp_types.h"
+#include "ortools/lp_data/lp_utils.h"
 #include "ortools/lp_data/scattered_vector.h"
 #include "ortools/lp_data/sparse.h"
 #include "ortools/util/bitset.h"
 #include "ortools/util/stats.h"
 
-#ifdef OMP
-#include <omp.h>
-#endif
-
-#include "ortools/lp_data/lp_utils.h"
-
 namespace operations_research {
 namespace glop {
 
@@ -275,8 +270,8 @@ void ReducedCosts::PerturbCosts() {
       case VariableType::UPPER_AND_LOWER_BOUNDED:
         // Here we don't necessarily maintain the dual-feasibility of a dual
         // feasible solution, however we can always shift the variable to its
-        // other bound (because it is boxed) to restore dual-feasiblity. This is
-        // done by MakeBoxedVariableDualFeasible() at the end of the dual
+        // other bound (because it is boxed) to restore dual-feasibility. This
+        // is done by MakeBoxedVariableDualFeasible() at the end of the dual
         // phase-I algorithm.
         if (objective > 0.0) {
           cost_perturbations_[col] = magnitude;
@@ -370,52 +365,30 @@ void ReducedCosts::ComputeReducedCosts() {
   }
   Fractional dual_residual_error(0.0);
   const ColIndex num_cols = matrix_.num_cols();
+  const ColIndex first_slack = num_cols - RowToColIndex(matrix_.num_rows());
 
   reduced_costs_.resize(num_cols, 0.0);
   const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow();
-#ifdef OMP
-  const int num_omp_threads = parameters_.num_omp_threads();
-#else
-  const int num_omp_threads = 1;
-#endif
-  if (num_omp_threads == 1) {
-    for (ColIndex col(0); col < num_cols; ++col) {
-      reduced_costs_[col] = objective_[col] + cost_perturbations_[col] -
-                            matrix_.ColumnScalarProduct(
-                                col, basic_objective_left_inverse_.values);
-
-      // We also compute the dual residual error y.B - c_B.
-      if (is_basic.IsSet(col)) {
-        dual_residual_error =
-            std::max(dual_residual_error, std::abs(reduced_costs_[col]));
-      }
-    }
-  } else {
-#ifdef OMP
-    // In the multi-threaded case, perform the same computation as in the
-    // single-threaded case above.
-    std::vector<Fractional> thread_local_dual_residual_error(num_omp_threads,
-                                                             0.0);
-    const int parallel_loop_size = num_cols.value();
-#pragma omp parallel for num_threads(num_omp_threads)
-    for (int i = 0; i < parallel_loop_size; i++) {
-      const ColIndex col(i);
-      reduced_costs_[col] = objective_[col] + objective_perturbation_[col] -
-                            matrix_.ColumnScalarProduct(
-                                col, basic_objective_left_inverse_.values);
-
-      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()],
-                     std::abs(reduced_costs_[col]));
-      }
+  for (ColIndex col(0); col < first_slack; ++col) {
+    reduced_costs_[col] =
+        objective_[col] + cost_perturbations_[col] -
+        matrix_.ColumnScalarProduct(col, basic_objective_left_inverse_.values);
+
+    // We also compute the dual residual error y.B - c_B.
+    if (is_basic.IsSet(col)) {
+      dual_residual_error =
+          std::max(dual_residual_error, std::abs(reduced_costs_[col]));
     }
-    // end of omp parallel for
-    for (int i = 0; i < num_omp_threads; i++) {
+  }
+  for (ColIndex col(first_slack); col < num_cols; ++col) {
+    reduced_costs_[col] = objective_[col] + cost_perturbations_[col] -
+                          basic_objective_left_inverse_[col - first_slack];
+
+    // We also compute the dual residual error y.B - c_B.
+    if (is_basic.IsSet(col)) {
       dual_residual_error =
-          std::max(dual_residual_error, thread_local_dual_residual_error[i]);
+          std::max(dual_residual_error, std::abs(reduced_costs_[col]));
     }
-#endif  // OMP
   }
 
   deterministic_time_ +=

ortools/glop/revised_simplex.cc

@@ -478,7 +478,7 @@ ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal(
       double min_distance = kInfinity;
       const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds();
       const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds();
-      double cost_delta = 0.0;
+      Fractional cost_delta = 0.0;
       for (ColIndex col(0); col < num_cols_; ++col) {
         cost_delta += solution_primal_ray_[col] * objective_[col];
         if (solution_primal_ray_[col] > 0 && upper_bounds[col] != kInfinity) {
@@ -586,10 +586,12 @@ ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal(
         // infeasibility lower than its corresponding residual error. Note that
         // we already adapt the tolerance like this during the simplex
         // execution.
-        const Fractional primal_tolerance = std::max(
-            primal_residual, parameters_.primal_feasibility_tolerance());
+        const Fractional primal_tolerance =
+            std::max(primal_residual,
+                     Fractional(parameters_.primal_feasibility_tolerance()));
         const Fractional dual_tolerance =
-            std::max(dual_residual, parameters_.dual_feasibility_tolerance());
+            std::max(dual_residual,
+                     Fractional(parameters_.dual_feasibility_tolerance()));
         const Fractional primal_infeasibility =
             variable_values_.ComputeMaximumPrimalInfeasibility();
         const Fractional dual_infeasibility =
@@ -2747,7 +2749,7 @@ Status RevisedSimplex::Polish(TimeLimit* time_limit) {
     const auto get_diff = [this](ColIndex col, Fractional old_value,
                                  Fractional new_value) {
       if (col >= integrality_scale_.size() || integrality_scale_[col] == 0.0) {
-        return 0.0;
+        return Fractional(0.0);
       }
       const Fractional s = integrality_scale_[col];
       return (std::abs(new_value * s - std::round(new_value * s)) -

ortools/glop/update_row.cc

@@ -304,14 +304,19 @@ void UpdateRow::ComputeUpdatesColumnWise() {
   non_zero_position_list_.resize(matrix_.num_cols().value());
   auto* non_zeros = non_zero_position_list_.data();
 
+  const ColIndex first_slack =
+      matrix_.num_cols() - RowToColIndex(matrix_.num_rows());
+
   const Fractional drop_tolerance = parameters_.drop_tolerance();
   const auto output_coeffs = coefficient_.view();
   const auto view = matrix_.view();
   const auto unit_row_left_inverse = unit_row_left_inverse_.values.const_view();
   for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) {
     // Coefficient of the column right inverse on the 'leaving_row'.
     const Fractional coeff =
-        view.ColumnScalarProduct(col, unit_row_left_inverse);
+        col >= first_slack
+            ? unit_row_left_inverse_[col - first_slack]
+            : view.ColumnScalarProduct(col, unit_row_left_inverse);
 
     // Nothing to do if 'coeff' is (almost) zero which does happen due to
     // sparsity. Note that it shouldn't be too bad to use a non-zero drop

ortools/glop/variable_values.cc

@@ -170,7 +170,7 @@ Fractional VariableValues::ComputeSumOfPrimalInfeasibilities() const {
   for (ColIndex col(0); col < num_cols; ++col) {
     const Fractional infeasibility =
         GetColInfeasibility(col, values, lower_bounds, upper_bounds);
-    sum += std::max(0.0, infeasibility);
+    sum += std::max(Fractional{0.0}, infeasibility);
   }
   return sum;
 }

ortools/lp_data/lp_parser.cc

@@ -235,6 +235,21 @@ bool LPParser::ParseConstraint(StringPiece constraint) {
   return true;
 }
 
+namespace {
+
+template <typename T>
+bool SimpleAtoFractional(absl::string_view str, T* value) {
+  if constexpr (std::is_same_v<T, double>) {
+    return absl::SimpleAtod(str, value);
+  } else if constexpr (std::is_same_v<T, float>) {
+    return absl::SimpleAtof(str, value);
+  } else {
+    static_assert(false, "Unsupported fractional type");
+    return false;
+  }
+}
+}  // namespace
+
 bool LPParser::SetVariableBounds(ColIndex col, Fractional lb, Fractional ub) {
   if (bounded_variables_.find(col) == bounded_variables_.end()) {
     // The variable was not bounded yet, thus reset its bounds.
@@ -250,7 +265,7 @@ bool LPParser::SetVariableBounds(ColIndex col, Fractional lb, Fractional ub) {
 }
 
 TokenType ConsumeToken(StringPiece* sp, std::string* consumed_name,
-                       double* consumed_coeff) {
+                       Fractional* consumed_coeff) {
   DCHECK(consumed_name != nullptr);
   DCHECK(consumed_coeff != nullptr);
   // We use LazyRE2 everywhere so that all the patterns are just compiled once
@@ -305,7 +320,7 @@ TokenType ConsumeToken(StringPiece* sp, std::string* consumed_name,
   static const LazyRE2 kValuePattern = {
       R"(\s*([0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?))"};
   if (RE2::Consume(sp, *kValuePattern, &coeff)) {
-    if (!absl::SimpleAtod(coeff, consumed_coeff)) {
+    if (!SimpleAtoFractional(coeff, consumed_coeff)) {
       // Note: If absl::SimpleAtod(), Consume(), and kValuePattern are correct,
       // this should never happen.
       LOG(ERROR) << "Text: " << coeff << " was matched by RE2 to be "

ortools/lp_data/lp_types.h

@@ -81,13 +81,13 @@ static inline double ToDouble(double f) { return f; }
 typedef double Fractional;
 
 // Range max for type Fractional. DBL_MAX for double for example.
-constexpr double kRangeMax = std::numeric_limits<double>::max();
+constexpr Fractional kRangeMax = std::numeric_limits<Fractional>::max();
 
 // Infinity for type Fractional.
-constexpr double kInfinity = std::numeric_limits<double>::infinity();
+constexpr Fractional kInfinity = std::numeric_limits<Fractional>::infinity();
 
 // Epsilon for type Fractional, i.e. the smallest e such that 1.0 + e != 1.0 .
-constexpr double kEpsilon = std::numeric_limits<double>::epsilon();
+constexpr Fractional kEpsilon = std::numeric_limits<Fractional>::epsilon();
 
 // Returns true if the given value is finite, that means for a double:
 // not a NaN and not +/- infinity.