Now bypass row/col scaling in Simplex NLA when only costs need to be unscaled

Author: jajhall

highs/lp_data/HighsLp.cpp

@@ -282,17 +282,24 @@ void HighsLp::unapplyScale() {
   // Unapply the scaling to the bounds, costs and matrix, and record
   // that it has been unapplied
   assert(scale.cost);
-  const double unscale_cost_mu = 1.0 / scale.cost;
-  for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) {
-    this->col_lower_[iCol] *= scale.col[iCol];
-    this->col_upper_[iCol] *= scale.col[iCol];
-    this->col_cost_[iCol] /= (unscale_cost_mu * scale.col[iCol]);
-  }
-  for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) {
-    this->row_lower_[iRow] /= scale.row[iRow];
-    this->row_upper_[iRow] /= scale.row[iRow];
+  const bool has_cost_scaling = scale.cost > 1.0;
+  const bool has_matrix_scaling = scale.col.size() && scale.row.size();
+  assert(has_cost_scaling || has_matrix_scaling);
+  if (has_matrix_scaling) {
+    for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) {
+      this->col_lower_[iCol] *= scale.col[iCol];
+      this->col_upper_[iCol] *= scale.col[iCol];
+      this->col_cost_[iCol] /= (scale.col[iCol] / scale.cost);
+    }
+    for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) {
+      this->row_lower_[iRow] /= scale.row[iRow];
+      this->row_upper_[iRow] /= scale.row[iRow];
+    }
+    this->a_matrix_.unapplyScale(scale);
+  } else {
+    for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) 
+      this->col_cost_[iCol] *= scale.cost;
   }
-  this->a_matrix_.unapplyScale(scale);
   this->is_scaled_ = false;
 }
 

highs/simplex/HSimplex.cpp

@@ -1056,18 +1056,31 @@ HighsStatus applyScalingToLpRow(HighsLp& lp, const HighsInt row,
 
 void simplexUnscaleSolution(HighsSolution& solution, const HighsScale& scale) {
   assert(scale.has_scaling);
-  assert(solution.col_value.size() == static_cast<size_t>(scale.num_col));
-  assert(solution.col_dual.size() == static_cast<size_t>(scale.num_col));
-  assert(solution.row_value.size() == static_cast<size_t>(scale.num_row));
-  assert(solution.row_dual.size() == static_cast<size_t>(scale.num_row));
-
-  for (HighsInt iCol = 0; iCol < scale.num_col; iCol++) {
-    solution.col_value[iCol] *= scale.col[iCol];
-    solution.col_dual[iCol] /= (scale.col[iCol] / scale.cost);
-  }
-  for (HighsInt iRow = 0; iRow < scale.num_row; iRow++) {
-    solution.row_value[iRow] /= scale.row[iRow];
-    solution.row_dual[iRow] *= (scale.row[iRow] * scale.cost);
+  const bool has_cost_scaling = scale.cost > 1.0;
+  const bool has_matrix_scaling = scale.col.size() && scale.row.size();
+  assert(has_cost_scaling || has_matrix_scaling);
+  HighsInt num_col = solution.col_value.size();
+  HighsInt num_row = solution.row_value.size();
+  if (has_matrix_scaling) {
+    assert(scale.num_col == num_col);
+    assert(scale.num_row == num_row);
+    assert(solution.col_value.size() == scale.col.size());
+    assert(solution.row_value.size() == scale.row.size());
+    assert(solution.col_dual.size() == scale.col.size());
+    assert(solution.row_dual.size() == scale.row.size());
+    for (HighsInt iCol = 0; iCol < num_col; iCol++) {
+      solution.col_value[iCol] *= scale.col[iCol];
+      solution.col_dual[iCol] /= (scale.col[iCol] / scale.cost);
+    }
+    for (HighsInt iRow = 0; iRow < num_row; iRow++) {
+      solution.row_value[iRow] /= scale.row[iRow];
+      solution.row_dual[iRow] *= (scale.row[iRow] * scale.cost);
+    }
+  } else {
+    for (HighsInt iCol = 0; iCol < num_col; iCol++) 
+      solution.col_dual[iCol] *= scale.cost;
+    for (HighsInt iRow = 0; iRow < num_row; iRow++) 
+      solution.row_dual[iRow] *= scale.cost;
   }
 }
 
@@ -1095,6 +1108,7 @@ void simplexScaleCost(const HighsOptions& options, HighsLp& lp) {
     cost_scale = pow(2.0, floor(log(cost_scale) / log(2.0) + 0.5));
     cost_scale = min(cost_scale, max_allowed_cost_scale);
   }
+  assert(cost_scale >= 1.0);
   if (cost_scale == 1.0) {
     highsLogDev(options.log_options, HighsLogType::kInfo,
 		"LP cost vector not scaled down: max cost is %g\n",