Merge branch 'latest' into fix-1869

FEATURES.md

@@ -20,6 +20,9 @@ Added minimal documentation of solvers and how simplex variants can be run
 
 Methods receiving matrix data where only small values are explicit zeros (so removed internally) are now silent and return HighsStatus::kOk (since internal matrix is exact)
 
+Now multiplying by pre-computed reciprocals rather than performing divisions in loops in simplex solver: LP performance improvement ~2.5%
+
+
 
 
 

src/simplex/HEkk.cpp

@@ -2126,6 +2126,7 @@ void HEkk::updateDualSteepestEdgeWeights(
 
   const double col_aq_scale = simplex_nla_.variableScaleFactor(variable_in);
   const double col_ap_scale = simplex_nla_.basicColScaleFactor(row_out);
+  const double inv_col_ap_scale = 1.0 / col_ap_scale;
 
   const bool DSE_check = false;
   HVector alt_dual_steepest_edge_column;
@@ -2165,7 +2166,7 @@ void HEkk::updateDualSteepestEdgeWeights(
     HighsInt DSE_array_count = 0;
     for (HighsInt iRow = 0; iRow < num_row; iRow++) {
       const double dual_steepest_edge_array_value =
-          dual_steepest_edge_array[iRow] / col_ap_scale;
+          dual_steepest_edge_array[iRow] * inv_col_ap_scale;
       if (dual_steepest_edge_array_value) DSE_array_count++;
       if (std::abs(dual_steepest_edge_array_value) >
               dse_column_value_tolerance ||
@@ -2194,7 +2195,7 @@ void HEkk::updateDualSteepestEdgeWeights(
              (int)alt_dual_steepest_edge_column.count, (int)DSE_array_count);
       for (HighsInt iRow = 0; iRow < num_row; iRow++) {
         const double dual_steepest_edge_array_value =
-            dual_steepest_edge_array[iRow] / col_ap_scale;
+            dual_steepest_edge_array[iRow] * inv_col_ap_scale;
         if (alt_dual_steepest_edge_column.array[iRow] != 0 &&
             dual_steepest_edge_array_value != 0) {
           const double dse_column_error =
@@ -2236,7 +2237,7 @@ void HEkk::updateDualSteepestEdgeWeights(
       double basic_col_scale = simplex_nla_.basicColScaleFactor(iRow);
       aa_iRow /= basic_col_scale;
       aa_iRow *= col_aq_scale;
-      dual_steepest_edge_array_value /= col_ap_scale;
+      dual_steepest_edge_array_value *= inv_col_ap_scale;
     }
     if (DSE_check) {
       const double pivotal_column_error =

src/simplex/HEkkDual.cpp

@@ -402,6 +402,7 @@ void HEkkDual::initialiseInstance() {
   solver_num_col = ekk_instance_.lp_.num_col_;
   solver_num_row = ekk_instance_.lp_.num_row_;
   solver_num_tot = solver_num_col + solver_num_row;
+  inv_solver_num_row = 1.0 / solver_num_row;
 
   a_matrix = &ekk_instance_.lp_.a_matrix_;
   simplex_nla = &ekk_instance_.simplex_nla_;
@@ -1276,7 +1277,7 @@ void HEkkDual::iterateTasks() {
   chooseRow();
 
   // Disable slice when too sparse
-  if (1.0 * row_ep.count / solver_num_row < 0.01) slice_PRICE = 0;
+  if (1.0 * row_ep.count * inv_solver_num_row < 0.01) slice_PRICE = 0;
 
   analysis->simplexTimerStart(Group1Clock);
   // #pragma omp parallel
@@ -1498,7 +1499,7 @@ void HEkkDual::chooseRow() {
   move_out = delta_primal < 0 ? -1 : 1;
   // Update the record of average row_ep (pi_p) density. This ignores
   // any BTRANs done for skipped candidates
-  const double local_row_ep_density = (double)row_ep.count / solver_num_row;
+  const double local_row_ep_density = (double)row_ep.count * inv_solver_num_row;
   ekk_instance_.updateOperationResultDensity(
       local_row_ep_density, ekk_instance_.info_.row_ep_density);
 }
@@ -1798,7 +1799,7 @@ void HEkkDual::chooseColumnSlice(HVector* row_ep) {
   analysis->simplexTimerStop(Chuzc0Clock);
 
   //  const HighsInt solver_num_row = ekk_instance_.lp_.num_row_;
-  const double local_density = 1.0 * row_ep->count / solver_num_row;
+  const double local_density = 1.0 * row_ep->count * inv_solver_num_row;
   bool use_col_price;
   bool use_row_price_w_switch;
   HighsSimplexInfo& info = ekk_instance_.info_;
@@ -1961,7 +1962,7 @@ void HEkkDual::updateFtran() {
                      analysis->pointer_serial_factor_clocks);
   if (analysis->analyse_simplex_summary_data)
     analysis->operationRecordAfter(kSimplexNlaFtran, col_aq);
-  const double local_col_aq_density = (double)col_aq.count / solver_num_row;
+  const double local_col_aq_density = (double)col_aq.count * inv_solver_num_row;
   ekk_instance_.updateOperationResultDensity(
       local_col_aq_density, ekk_instance_.info_.col_aq_density);
   // Save the pivot value computed column-wise - used for numerical checking
@@ -2002,7 +2003,8 @@ void HEkkDual::updateFtranBFRT() {
   if (time_updateFtranBFRT) {
     analysis->simplexTimerStop(FtranBfrtClock);
   }
-  const double local_col_BFRT_density = (double)col_BFRT.count / solver_num_row;
+  const double local_col_BFRT_density =
+      (double)col_BFRT.count * inv_solver_num_row;
   ekk_instance_.updateOperationResultDensity(
       local_col_BFRT_density, ekk_instance_.info_.col_BFRT_density);
 }
@@ -2043,7 +2045,7 @@ void HEkkDual::updateFtranDSE(HVector* DSE_Vector) {
     analysis->operationRecordAfter(kSimplexNlaFtranDse, *DSE_Vector);
   analysis->simplexTimerStop(FtranDseClock);
   const double local_row_DSE_density =
-      (double)DSE_Vector->count / solver_num_row;
+      (double)DSE_Vector->count * inv_solver_num_row;
   ekk_instance_.updateOperationResultDensity(
       local_row_DSE_density, ekk_instance_.info_.row_DSE_density);
 }

src/simplex/HEkkDual.h

@@ -388,6 +388,7 @@ class HEkkDual {
   HighsInt solver_num_row;
   HighsInt solver_num_col;
   HighsInt solver_num_tot;
+  double inv_solver_num_row;  // 1.0 / solver_num_row
 
   const HighsSparseMatrix* a_matrix;
   const HSimplexNla* simplex_nla;

src/util/HFactor.cpp

@@ -211,6 +211,7 @@ void HFactor::setupGeneral(
   num_row = num_row_;
   num_col = num_col_;
   num_basic = num_basic_;
+  inv_num_row = 1.0 / num_row;
   this->a_matrix_valid = true;
   a_start = a_start_;
   a_index = a_index_;
@@ -1543,7 +1544,7 @@ void HFactor::ftranL(HVector& rhs, const double expected_density,
   }
 
   // Determine style of solve
-  double current_density = 1.0 * rhs.count / num_row;
+  double current_density = 1.0 * rhs.count * inv_num_row;
   const bool sparse_solve = rhs.count < 0 || current_density > kHyperCancel ||
                             expected_density > kHyperFtranL;
   if (sparse_solve) {
@@ -1591,7 +1592,7 @@ void HFactor::btranL(HVector& rhs, const double expected_density,
   factor_timer.start(FactorBtranLower, factor_timer_clock_pointer);
 
   // Determine style of solve
-  const double current_density = 1.0 * rhs.count / num_row;
+  const double current_density = 1.0 * rhs.count * inv_num_row;
   const bool sparse_solve = rhs.count < 0 || current_density > kHyperCancel ||
                             expected_density > kHyperBtranL;
   if (sparse_solve) {
@@ -1666,7 +1667,7 @@ void HFactor::ftranU(HVector& rhs, const double expected_density,
   // The regular part
   //
   // Determine style of solve
-  const double current_density = 1.0 * rhs.count / num_row;
+  const double current_density = 1.0 * rhs.count * inv_num_row;
   const bool sparse_solve = rhs.count < 0 || current_density > kHyperCancel ||
                             expected_density > kHyperFtranU;
   if (sparse_solve) {
@@ -1720,7 +1721,7 @@ void HFactor::ftranU(HVector& rhs, const double expected_density,
         rhs_synthetic_tick * 15 + (u_pivot_count - num_row) * 10;
     factor_timer.stop(use_clock, factor_timer_clock_pointer);
     if (report_ftran_upper_sparse) {
-      const double final_density = 1.0 * rhs.count / num_row;
+      const double final_density = 1.0 * rhs.count * inv_num_row;
       printf(
           "FactorFtranUpperSps: expected_density = %10.4g; current_density = "
           "%10.4g; final_density = %10.4g\n",
@@ -1773,7 +1774,7 @@ void HFactor::btranU(HVector& rhs, const double expected_density,
   // The regular part
   //
   // Determine style of solve
-  const double current_density = 1.0 * rhs.count / num_row;
+  const double current_density = 1.0 * rhs.count * inv_num_row;
   const bool sparse_solve = rhs.count < 0 || current_density > kHyperCancel ||
                             expected_density > kHyperBtranU;
   if (sparse_solve) {

src/util/HFactor.h

@@ -347,6 +347,7 @@ class HFactor {
   HighsInt num_row;
   HighsInt num_col;
   HighsInt num_basic;
+  double inv_num_row;  // 1.0/num_row
 
  private:
   bool a_matrix_valid;

src/util/HighsSparseMatrix.cpp

@@ -1383,6 +1383,7 @@ void HighsSparseMatrix::priceByRowWithSwitch(
   assert(HighsInt(result.size) == this->num_col_);
   assert(HighsInt(result.index.size()) == this->num_col_);
   if (expected_density <= kHyperPriceDensity) {
+    double inv_num_col = 1.0 / this->num_col_;
     for (HighsInt ix = next_index; ix < column.count; ix++) {
       HighsInt iRow = column.index[ix];
       // Determine whether p_end_ or the next start_ ends the loop
@@ -1394,7 +1395,7 @@ void HighsSparseMatrix::priceByRowWithSwitch(
       }
       // Possibly switch to standard row-wise price
       HighsInt row_num_nz = to_iEl - this->start_[iRow];
-      double local_density = (1.0 * result.count) / this->num_col_;
+      double local_density = (1.0 * result.count) * inv_num_col;
       bool switch_to_dense = result.count + row_num_nz >= this->num_col_ ||
                              local_density > switch_density;
       if (switch_to_dense) break;