Fix 2087

check/TestIpm.cpp

@@ -134,3 +134,22 @@ TEST_CASE("test-1966", "[highs_ipm]") {
     printf("Sum   dual infeasibilities = %g\n", info.sum_dual_infeasibilities);
   }
 }
+
+TEST_CASE("test-2087", "[highs_ipm]") {
+  // Make sure that presolve is performed when re-solving using IPM,
+  // since optimal basis cannot be used, and ensure that the offset is
+  // used in IPX
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  // Use shell since it yields an offset after presolve
+  std::string model = "shell.mps";
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/" + model;
+  h.readModel(filename);
+
+  h.setOptionValue("solver", kIpmString);
+  h.run();
+  const HighsInt first_ipm_iteration_count = h.getInfo().ipm_iteration_count;
+
+  h.run();
+  REQUIRE(first_ipm_iteration_count == h.getInfo().ipm_iteration_count);
+}

check/TestIpx.cpp

@@ -22,6 +22,7 @@ using Int = ipxint;
 
 constexpr HighsInt num_var = 12;
 constexpr HighsInt num_constr = 9;
+const double offset = 0;
 const double obj[] = {-0.2194, 0.0, 0.0,   0.0,     0.0, 0.0,
                       0.0,     0.0, -0.32, -0.5564, 0.6, -0.48};
 const double lb[num_var] = {0.0};
@@ -47,8 +48,8 @@ TEST_CASE("test-ipx", "[highs_ipx]") {
   lps.SetParameters(parameters);
 
   // Solve the LP.
-  Int load_status = lps.LoadModel(num_var, obj, lb, ub, num_constr, Ap, Ai, Ax,
-                                  rhs, constr_type);
+  Int load_status = lps.LoadModel(num_var, offset, obj, lb, ub, num_constr, Ap,
+                                  Ai, Ax, rhs, constr_type);
   REQUIRE(load_status == 0);
 
   highs::parallel::initialize_scheduler();

src/ipm/IpxWrapper.cpp

@@ -150,18 +150,19 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   lps.SetCallback(&callback);
 
   ipx::Int num_col, num_row;
+  double offset;
   std::vector<ipx::Int> Ap, Ai;
   std::vector<double> objective, col_lb, col_ub, Av, rhs;
   std::vector<char> constraint_type;
-  fillInIpxData(lp, num_col, num_row, objective, col_lb, col_ub, Ap, Ai, Av,
-                rhs, constraint_type);
+  fillInIpxData(lp, num_col, num_row, offset, objective, col_lb, col_ub, Ap, Ai,
+                Av, rhs, constraint_type);
   highsLogUser(options.log_options, HighsLogType::kInfo,
                "IPX model has %" HIGHSINT_FORMAT " rows, %" HIGHSINT_FORMAT
                " columns and %" HIGHSINT_FORMAT " nonzeros\n",
                num_row, num_col, Ap[num_col]);
 
   ipx::Int load_status = lps.LoadModel(
-      num_col, objective.data(), col_lb.data(), col_ub.data(), num_row,
+      num_col, offset, objective.data(), col_lb.data(), col_ub.data(), num_row,
       Ap.data(), Ai.data(), Av.data(), rhs.data(), constraint_type.data());
 
   if (load_status) {
@@ -387,10 +388,10 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
 }
 
 void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row,
-                   std::vector<double>& obj, std::vector<double>& col_lb,
-                   std::vector<double>& col_ub, std::vector<ipx::Int>& Ap,
-                   std::vector<ipx::Int>& Ai, std::vector<double>& Ax,
-                   std::vector<double>& rhs,
+                   double& offset, std::vector<double>& obj,
+                   std::vector<double>& col_lb, std::vector<double>& col_ub,
+                   std::vector<ipx::Int>& Ap, std::vector<ipx::Int>& Ai,
+                   std::vector<double>& Ax, std::vector<double>& rhs,
                    std::vector<char>& constraint_type) {
   num_col = lp.num_col_;
   num_row = lp.num_row_;
@@ -517,6 +518,7 @@ void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row,
     col_ub[lp.num_col_ + slack] = lp.row_upper_[row];
   }
 
+  offset = HighsInt(lp.sense_) * lp.offset_;
   obj.resize(num_col);
   for (HighsInt col = 0; col < lp.num_col_; col++) {
     obj[col] = (HighsInt)lp.sense_ * lp.col_cost_[col];

src/ipm/IpxWrapper.h

@@ -28,10 +28,10 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
                        HighsCallback& callback);
 
 void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row,
-                   std::vector<double>& obj, std::vector<double>& col_lb,
-                   std::vector<double>& col_ub, std::vector<ipx::Int>& Ap,
-                   std::vector<ipx::Int>& Ai, std::vector<double>& Ax,
-                   std::vector<double>& rhs,
+                   double& offset, std::vector<double>& obj,
+                   std::vector<double>& col_lb, std::vector<double>& col_ub,
+                   std::vector<ipx::Int>& Ap, std::vector<ipx::Int>& Ai,
+                   std::vector<double>& Ax, std::vector<double>& rhs,
                    std::vector<char>& constraint_type);
 
 HighsStatus reportIpxSolveStatus(const HighsOptions& options,

src/ipm/ipx/ipx_c.cc

@@ -9,7 +9,8 @@ ipxint ipx_load_model(void* self, ipxint num_var, const double* obj,
                       const ipxint* Ap, const ipxint* Ai, const double* Ax,
                       const double* rhs, const char* constr_type) {
     LpSolver* solver = static_cast<LpSolver*>(self);
-    return solver->LoadModel(num_var, obj, lb, ub, num_constr, Ap, Ai, Ax, rhs,
+    const double offset = 0;
+    return solver->LoadModel(num_var, offset, obj, lb, ub, num_constr, Ap, Ai, Ax, rhs,
                              constr_type);
 }
 

src/ipm/ipx/iterate.cc

@@ -617,8 +617,8 @@ void Iterate::ComputeObjectives() const {
     if (postprocessed_) {
         // Compute objective values as defined for the LP model.
         offset_ = 0.0;
-        pobjective_ = Dot(c, x_);
-        dobjective_ = Dot(b, y_);
+        pobjective_ = model_.offset() + Dot(c, x_);
+        dobjective_ = model_.offset() + Dot(b, y_);
         for (Int j = 0; j < n+m; j++) {
             if (std::isfinite(lb[j]))
                 dobjective_ += lb[j] * zl_[j];
@@ -630,7 +630,7 @@ void Iterate::ComputeObjectives() const {
         // (after fixing and implying variables). The offset is such that
         // pobjective_ + offset_ is the primal objective after postprocessing.
         offset_ = 0.0;
-        pobjective_ = 0.0;
+        pobjective_ = model_.offset();
         for (Int j = 0; j < n+m; j++) {
             if (StateOf(j) != State::fixed)
                 pobjective_ += c[j] * x_[j];
@@ -643,7 +643,7 @@ void Iterate::ComputeObjectives() const {
                 offset_ += (zl_[j]-zu_[j]) * x_[j];
             }
         }
-        dobjective_ = Dot(b, y_);
+        dobjective_ = model_.offset() + Dot(b, y_);
         for (Int j = 0; j < n+m; j++) {
             if (has_barrier_lb(j))
                 dobjective_ += lb[j] * zl_[j];

src/ipm/ipx/lp_solver.cc

@@ -12,13 +12,14 @@
 
 namespace ipx {
 
-Int LpSolver::LoadModel(Int num_var, const double* obj, const double* lb,
+Int LpSolver::LoadModel(Int num_var, const double offset,
+			const double* obj, const double* lb,
                         const double* ub, Int num_constr, const Int* Ap,
                         const Int* Ai, const double* Ax, const double* rhs,
                         const char* constr_type) {
     ClearModel();
     Int errflag = model_.Load(control_, num_constr, num_var, Ap, Ai, Ax, rhs,
-                              constr_type, obj, lb, ub);
+                              constr_type, offset, obj, lb, ub);
     model_.GetInfo(&info_);
     return errflag;
 }

src/ipm/ipx/lp_solver.h

@@ -28,7 +28,8 @@ class LpSolver {
     //  IPX_ERROR_invalid_dimension
     //  IPX_ERROR_invalid_matrix
     //  IPX_ERROR_invalid_vector
-    Int LoadModel(Int num_var, const double* obj, const double* lb,
+    Int LoadModel(Int num_var, const double offset,
+		  const double* obj, const double* lb,
                   const double* ub, Int num_constr, const Int* Ap,
                   const Int* Ai, const double* Ax, const double* rhs,
                   const char* constr_type);

src/ipm/ipx/model.cc

@@ -9,11 +9,11 @@ namespace ipx {
 
 Int Model::Load(const Control& control, Int num_constr, Int num_var,
                 const Int* Ap, const Int* Ai, const double* Ax,
-                const double* rhs, const char* constr_type, const double* obj,
-                const double* lbuser, const double* ubuser) {
+                const double* rhs, const char* constr_type, const double offset,
+		const double* obj, const double* lbuser, const double* ubuser) {
     clear();
     Int errflag = CopyInput(num_constr, num_var, Ap, Ai, Ax, rhs, constr_type,
-                            obj, lbuser, ubuser);
+                            offset, obj, lbuser, ubuser);
     if (errflag)
         return errflag;
     std::stringstream h_logging_stream;
@@ -336,8 +336,8 @@ void Model::EvaluateInteriorSolution(const Vector& x_solver,
     presidual = std::max(presidual, Infnorm(ru));
     double dresidual = Infnorm(rc);
 
-    double pobjective = Dot(scaled_obj_, x);
-    double dobjective = Dot(scaled_rhs_, y);
+    double pobjective = offset_ + Dot(scaled_obj_, x);
+    double dobjective = offset_ + Dot(scaled_rhs_, y);
     for (Int j = 0; j < num_var_; j++) {
         if (std::isfinite(scaled_lbuser_[j]))
             dobjective += scaled_lbuser_[j] * zl[j];
@@ -541,7 +541,7 @@ static Int CheckMatrix(Int m, Int n, const Int *Ap, const Int *Ai, const double
 
 Int Model::CopyInput(Int num_constr, Int num_var, const Int* Ap, const Int* Ai,
                      const double* Ax, const double* rhs,
-                     const char* constr_type, const double* obj,
+                     const char* constr_type,  const double offset, const double* obj,
                      const double* lbuser, const double* ubuser) {
     if (!(Ap && Ai && Ax && rhs && constr_type && obj && lbuser && ubuser)) {
         return IPX_ERROR_argument_null;
@@ -569,6 +569,7 @@ Int Model::CopyInput(Int num_constr, Int num_var, const Int* Ap, const Int* Ai,
             boxed_vars_.push_back(j);
     }
     constr_type_ = std::vector<char>(constr_type, constr_type+num_constr);
+    offset_ = offset;
     scaled_obj_ = Vector(obj, num_var);
     scaled_rhs_ = Vector(rhs, num_constr);
     scaled_lbuser_ = Vector(lbuser, num_var);

src/ipm/ipx/model.h

@@ -58,8 +58,8 @@ class Model {
     //  IPX_ERROR_invalid_vector
     Int Load(const Control& control, Int num_constr, Int num_var,
              const Int* Ap, const Int* Ai, const double* Ax,
-             const double* rhs, const char* constr_type, const double* obj,
-             const double* lbuser, const double* ubuser);
+             const double* rhs, const char* constr_type, const double offset,
+	     const double* obj, const double* lbuser, const double* ubuser);
     // Performs Flippo's test for deciding dualization
     bool filippoDualizationTest() const;
     // Writes statistics of input data and preprocessing to @info.
@@ -93,6 +93,9 @@ class Model {
     const SparseMatrix& AI() const { return AI_; }
     const SparseMatrix& AIt() const { return AIt_; }
 
+    // Returns the offset
+    const double offset() const { return offset_; }
+
     // Returns a reference to a model vector.
     const Vector& b() const { return b_; }
     const Vector& c() const { return c_; }
@@ -207,7 +210,7 @@ class Model {
     //  IPX_ERROR_invalid_vector
     Int CopyInput(Int num_constr, Int num_var, const Int* Ap, const Int* Ai,
                   const double* Ax, const double* rhs, const char* constr_type,
-                  const double* obj, const double* lbuser,
+                  const double offset, const double* obj, const double* lbuser,
                   const double* ubuser);
 
     // Scales A_, scaled_obj_, scaled_rhs_, scaled_lbuser_ and scaled_ubuser_
@@ -376,6 +379,7 @@ class Model {
     std::vector<char> constr_type_;
     double norm_obj_{0.0};        // Infnorm(obj) as given by user
     double norm_rhs_{0.0};        // Infnorm(rhs,lb,ub) as given by user
+    double offset_;
     Vector scaled_obj_;
     Vector scaled_rhs_;
     Vector scaled_lbuser_;

src/lp_data/Highs.cpp

@@ -1260,10 +1260,16 @@ HighsStatus Highs::solve() {
 
   const bool unconstrained_lp = incumbent_lp.a_matrix_.numNz() == 0;
   assert(incumbent_lp.num_row_ || unconstrained_lp);
-  if (basis_.valid || options_.presolve == kHighsOffString ||
-      unconstrained_lp) {
+  // Even if options_.solver == kHighsChooseString in isolation will,
+  // untimately lead to a choice between simplex and IPM, if a basis
+  // is available, simplex should surely be chosen.
+  const bool solver_will_use_basis = options_.solver == kSimplexString ||
+                                     options_.solver == kHighsChooseString;
+  if ((basis_.valid || options_.presolve == kHighsOffString ||
+       unconstrained_lp) &&
+      solver_will_use_basis) {
     // There is a valid basis for the problem, presolve is off, or LP
-    // has no constraint matrix
+    // has no constraint matrix, and the solver will use the basis
     ekk_instance_.lp_name_ =
         "LP without presolve, or with basis, or unconstrained";
     // If there is a valid HiGHS basis, refine any status values that

src/presolve/ICrashX.cpp

@@ -20,12 +20,13 @@ HighsStatus callCrossover(const HighsOptions& options, const HighsLp& lp,
                           HighsModelStatus& model_status, HighsInfo& highs_info,
                           HighsCallback& highs_callback) {
   ipx::Int num_col, num_row;
+  double offset;
   std::vector<ipx::Int> Ap, Ai;
   std::vector<double> objective, col_lb, col_ub, Av, rhs;
   std::vector<char> constraint_type;
 
-  fillInIpxData(lp, num_col, num_row, objective, col_lb, col_ub, Ap, Ai, Av,
-                rhs, constraint_type);
+  fillInIpxData(lp, num_col, num_row, offset, objective, col_lb, col_ub, Ap, Ai,
+                Av, rhs, constraint_type);
   // if (res != IpxStatus::OK) return HighsStatus::kError;
 
   const HighsLogOptions& log_options = options.log_options;
@@ -55,7 +56,7 @@ HighsStatus callCrossover(const HighsOptions& options, const HighsLp& lp,
   lps.SetCallback(&highs_callback);
 
   ipx::Int load_status = lps.LoadModel(
-      num_col, objective.data(), col_lb.data(), col_ub.data(), num_row,
+      num_col, offset, objective.data(), col_lb.data(), col_ub.data(), num_row,
       Ap.data(), Ai.data(), Av.data(), rhs.data(), constraint_type.data());
   if (load_status != 0) {
     highsLogUser(log_options, HighsLogType::kError,