diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index cb2d1a3419..dd59b96aea 100644 --- a/check/TestIpm.cpp +++ b/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); +} diff --git a/check/TestIpx.cpp b/check/TestIpx.cpp index 4866c746d3..a53a199710 100644 --- a/check/TestIpx.cpp +++ b/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(); diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index ccd7e4b3dd..7a33b68199 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/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 Ap, Ai; std::vector objective, col_lb, col_ub, Av, rhs; std::vector 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& obj, std::vector& col_lb, - std::vector& col_ub, std::vector& Ap, - std::vector& Ai, std::vector& Ax, - std::vector& rhs, + double& offset, std::vector& obj, + std::vector& col_lb, std::vector& col_ub, + std::vector& Ap, std::vector& Ai, + std::vector& Ax, std::vector& rhs, std::vector& 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]; diff --git a/src/ipm/IpxWrapper.h b/src/ipm/IpxWrapper.h index 28a46fe202..99486c2107 100644 --- a/src/ipm/IpxWrapper.h +++ b/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& obj, std::vector& col_lb, - std::vector& col_ub, std::vector& Ap, - std::vector& Ai, std::vector& Ax, - std::vector& rhs, + double& offset, std::vector& obj, + std::vector& col_lb, std::vector& col_ub, + std::vector& Ap, std::vector& Ai, + std::vector& Ax, std::vector& rhs, std::vector& constraint_type); HighsStatus reportIpxSolveStatus(const HighsOptions& options, diff --git a/src/ipm/ipx/ipx_c.cc b/src/ipm/ipx/ipx_c.cc index e70f1dbfed..33c2817f99 100644 --- a/src/ipm/ipx/ipx_c.cc +++ b/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(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); } diff --git a/src/ipm/ipx/iterate.cc b/src/ipm/ipx/iterate.cc index c407eb3de4..030620ea0d 100644 --- a/src/ipm/ipx/iterate.cc +++ b/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]; diff --git a/src/ipm/ipx/lp_solver.cc b/src/ipm/ipx/lp_solver.cc index 145930c75b..2276ec4e62 100644 --- a/src/ipm/ipx/lp_solver.cc +++ b/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; } diff --git a/src/ipm/ipx/lp_solver.h b/src/ipm/ipx/lp_solver.h index 18eb70606f..a7a2da4bd1 100644 --- a/src/ipm/ipx/lp_solver.h +++ b/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); diff --git a/src/ipm/ipx/model.cc b/src/ipm/ipx/model.cc index 94d251c2e4..321a419d82 100644 --- a/src/ipm/ipx/model.cc +++ b/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(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); diff --git a/src/ipm/ipx/model.h b/src/ipm/ipx/model.h index e57fe4bd0e..a06dec3ec0 100644 --- a/src/ipm/ipx/model.h +++ b/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 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_; diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index 1006ce43a0..61320a1b4a 100644 --- a/src/lp_data/Highs.cpp +++ b/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 diff --git a/src/presolve/ICrashX.cpp b/src/presolve/ICrashX.cpp index b4eab7d1e3..45e0c9cc5e 100644 --- a/src/presolve/ICrashX.cpp +++ b/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 Ap, Ai; std::vector objective, col_lb, col_ub, Av, rhs; std::vector 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,