diff --git a/check/TestIpm.cpp b/check/TestIpm.cpp index 911d01a67a..cb2d1a3419 100644 --- a/check/TestIpm.cpp +++ b/check/TestIpm.cpp @@ -81,3 +81,56 @@ TEST_CASE("test-analytic-centre-box", "[highs_ipm]") { REQUIRE(solution_norm < 1e-6); if (dev_run) printf("Analytic centre solution norm is %g\n", solution_norm); } + +TEST_CASE("test-1966", "[highs_ipm]") { + Highs highs; + highs.setOptionValue("output_flag", dev_run); + const HighsInfo& info = highs.getInfo(); + HighsLp lp; + lp.num_col_ = 2; + lp.num_row_ = 2; + lp.col_cost_ = {2, -1}; + lp.col_lower_ = {0, 0}; + lp.col_upper_ = {kHighsInf, kHighsInf}; + lp.row_lower_ = {-kHighsInf, 2}; + lp.row_upper_ = {1, kHighsInf}; + lp.a_matrix_.start_ = {0, 2, 4}; + lp.a_matrix_.index_ = {0, 1, 0, 1}; + lp.a_matrix_.value_ = {1, -1, 1, -1}; + highs.passModel(lp); + highs.setOptionValue("solver", kIpmString); + highs.setOptionValue("presolve", kHighsOffString); + + if (dev_run) printf("\nWith default residual tolerances\n"); + highs.run(); + if (dev_run) { + highs.writeSolution("", kSolutionStylePretty); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } + highs.clearSolver(); + + if (dev_run) printf("\nWith infinite residual tolerances\n"); + highs.setOptionValue("primal_residual_tolerance", 1e30); + highs.setOptionValue("dual_residual_tolerance", 1e30); + highs.run(); + if (dev_run) { + highs.writeSolution("", kSolutionStylePretty); + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } +} diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp index 5202a1dee5..572a478cf2 100644 --- a/check/TestLpSolvers.cpp +++ b/check/TestLpSolvers.cpp @@ -438,3 +438,35 @@ TEST_CASE("dual-objective-upper-bound", "[highs_lp_solver]") { if (dev_run) printf("\nOptimal objective value error = %g\n", error); REQUIRE(error < 1e-10); } + +TEST_CASE("blending-lp-ipm", "[highs_lp_solver]") { + Highs highs; + highs.setOptionValue("output_flag", dev_run); + HighsLp lp; + lp.num_col_ = 2; + lp.num_row_ = 2; + lp.col_cost_ = {-8, -10}; + lp.col_lower_ = {0, 0}; + lp.col_upper_ = {kHighsInf, kHighsInf}; + lp.row_lower_ = {-kHighsInf, -kHighsInf}; + lp.row_upper_ = {80, 120}; + lp.a_matrix_.start_ = {0, 2, 4}; + lp.a_matrix_.index_ = {0, 1, 0, 1}; + lp.a_matrix_.value_ = {1, 1, 2, 4}; + highs.passModel(lp); + highs.setOptionValue("solver", kIpmString); + highs.setOptionValue("presolve", kHighsOffString); + highs.run(); + HighsInfo info = highs.getInfo(); + if (dev_run) { + printf("Num primal infeasibilities = %d\n", + int(info.num_primal_infeasibilities)); + printf("Max primal infeasibility = %g\n", info.max_primal_infeasibility); + printf("Sum primal infeasibilities = %g\n", + info.sum_primal_infeasibilities); + printf("Num dual infeasibilities = %d\n", + int(info.num_dual_infeasibilities)); + printf("Max dual infeasibility = %g\n", info.max_dual_infeasibility); + printf("Sum dual infeasibilities = %g\n", info.sum_dual_infeasibilities); + } +} diff --git a/check/TestNames.cpp b/check/TestNames.cpp index f0d11fc4a8..4aceee0ca6 100644 --- a/check/TestNames.cpp +++ b/check/TestNames.cpp @@ -120,13 +120,12 @@ TEST_CASE("highs-names", "[highs_names]") { } TEST_CASE("highs-model-name", "[model_names]") { - Highs highs; const HighsLp& lp = highs.getLp(); std::string name = lp.model_name_; REQUIRE(name == ""); - + highs.passModelName("new_name"); name = lp.model_name_; REQUIRE(name == "new_name"); diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp index c420839a5d..172df328c4 100644 --- a/src/lp_data/Highs.cpp +++ b/src/lp_data/Highs.cpp @@ -1959,10 +1959,6 @@ HighsStatus Highs::getReducedColumn(const HighsInt col, double* col_vector, HighsStatus Highs::getKappa(double& kappa, const bool exact, const bool report) { - printf( - "Highs::getKappa basis_.valid = %d, ekk_instance_.status_.has_invert = " - "%d\n", - int(basis_.valid), int(ekk_instance_.status_.has_invert)); if (!ekk_instance_.status_.has_invert) return invertRequirementError("getBasisInverseRow"); kappa = ekk_instance_.computeBasisCondition(this->model_.lp_, exact, report); diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 208bcd1802..bf853a4839 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -389,6 +389,8 @@ struct HighsOptionsStruct { bool less_infeasible_DSE_choose_row; bool use_original_HFactor_logic; bool run_centring; + double primal_residual_tolerance; + double dual_residual_tolerance; HighsInt max_centring_steps; double centring_ratio_tolerance; @@ -524,6 +526,8 @@ struct HighsOptionsStruct { run_centring(false), max_centring_steps(0), centring_ratio_tolerance(0.0), + primal_residual_tolerance(0.0), + dual_residual_tolerance(0.0), icrash(false), icrash_dualize(false), icrash_strategy(""), @@ -1378,6 +1382,16 @@ class HighsOptions : public HighsOptionsStruct { advanced, ¢ring_ratio_tolerance, 0, 100, kHighsInf); records.push_back(record_double); + record_double = new OptionRecordDouble( + "primal_residual_tolerance", "Primal residual tolerance", advanced, + &primal_residual_tolerance, 1e-10, 1e-7, kHighsInf); + records.push_back(record_double); + + record_double = new OptionRecordDouble( + "dual_residual_tolerance", "Dual residual tolerance", advanced, + &dual_residual_tolerance, 1e-10, 1e-7, kHighsInf); + records.push_back(record_double); + // Set up the log_options aliases log_options.clear(); log_options.log_stream = diff --git a/src/lp_data/HighsSolution.cpp b/src/lp_data/HighsSolution.cpp index 29e40168a0..27cd16153e 100644 --- a/src/lp_data/HighsSolution.cpp +++ b/src/lp_data/HighsSolution.cpp @@ -282,7 +282,8 @@ void getKktFailures(const HighsOptions& options, const HighsLp& lp, relative_primal_infeasibility, dual_infeasibility, value_residual); if (!status_value_ok) highsLogUser(options.log_options, HighsLogType::kError, - "getKktFailures: %s %d status-value error: [%g; %g; %g] has " + "getKktFailures: %s %d status-value error: [%23.18g; " + "%23.18g; %23.18g] has " "residual %g\n", iVar < lp.num_col_ ? "Column" : "Row ", iVar < lp.num_col_ ? int(iVar) : int(iVar - lp.num_col_), @@ -531,14 +532,20 @@ bool getVariableKktFailures(const double primal_feasibility_tolerance, // If the variable is basic, then consider it not to be at a bound // so that any dual value yields an infeasibility value if (*status_pointer == HighsBasisStatus::kBasic) at_a_bound = false; + // With very large values, accuracy is lost in adding/subtracting + // the feasibility tolerance from the bounds, so skip if this may + // occur + // // Check that kLower and kUpper are consistent with value and // bounds - for debugging QP basis errors if (*status_pointer == HighsBasisStatus::kLower) { - status_value_ok = value >= lower - primal_feasibility_tolerance && - value <= lower + primal_feasibility_tolerance; + if (std::fabs(lower) / primal_feasibility_tolerance < 1e-16) + status_value_ok = value >= lower - primal_feasibility_tolerance && + value <= lower + primal_feasibility_tolerance; } else if (*status_pointer == HighsBasisStatus::kUpper) { - status_value_ok = value >= upper - primal_feasibility_tolerance && - value <= upper + primal_feasibility_tolerance; + if (std::fabs(upper) / primal_feasibility_tolerance < 1e-16) + status_value_ok = value >= upper - primal_feasibility_tolerance && + value <= upper + primal_feasibility_tolerance; } } if (at_a_bound) { @@ -1369,6 +1376,73 @@ void accommodateAlienBasis(HighsLpSolverObject& solver_object) { assert(num_basic_variables == num_row); } +void correctResiduals(HighsLpSolverObject& solver_object) { + HighsOptions& options = solver_object.options_; + HighsSolution& solution = solver_object.solution_; + const bool& have_primal_solution = solution.value_valid; + const bool& have_dual_solution = solution.dual_valid; + assert(have_primal_solution); + assert(have_dual_solution); + HighsLp& lp = solver_object.lp_; + std::vector check_row_value; + std::vector check_col_dual; + lp.a_matrix_.productQuad(check_row_value, solution.col_value); + if (have_dual_solution) { + lp.a_matrix_.productTransposeQuad(check_col_dual, solution.row_dual); + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) + check_col_dual[iCol] -= lp.col_cost_[iCol]; + } + double norm_primal_residual = 0; + HighsInt num_primal_correction = 0; + double max_primal_correction = 0; + double sum_primal_correction = 0; + double use_primal_residual_tolerance = options.primal_residual_tolerance; + + double norm_dual_residual = 0; + HighsInt num_dual_correction = 0; + double max_dual_correction = 0; + double sum_dual_correction = 0; + double use_dual_residual_tolerance = options.dual_residual_tolerance; + + for (HighsInt iRow = 0; iRow < lp.num_row_; iRow++) { + double primal_residual = check_row_value[iRow] - solution.row_value[iRow]; + double abs_primal_residual = std::fabs(primal_residual); + if (abs_primal_residual > use_primal_residual_tolerance) { + solution.row_value[iRow] += primal_residual; + double check_primal_residual = + std::fabs(check_row_value[iRow] - solution.row_value[iRow]); + assert(check_primal_residual < 1e-10); + num_primal_correction++; + max_primal_correction = + std::max(abs_primal_residual, max_primal_correction); + sum_primal_correction += abs_primal_residual; + } + norm_primal_residual = std::max(abs_primal_residual, norm_primal_residual); + } + if (have_dual_solution) { + for (HighsInt iCol = 0; iCol < lp.num_col_; iCol++) { + double dual_residual = check_col_dual[iCol] + solution.col_dual[iCol]; + double abs_dual_residual = std::fabs(dual_residual); + if (abs_dual_residual > use_dual_residual_tolerance) { + solution.col_dual[iCol] -= dual_residual; + num_dual_correction++; + max_dual_correction = std::max(abs_dual_residual, max_dual_correction); + sum_dual_correction += abs_dual_residual; + } + norm_dual_residual = std::max(abs_dual_residual, norm_dual_residual); + } + } + + if (num_primal_correction > 0 || num_dual_correction > 0) + highsLogUser( + options.log_options, HighsLogType::kWarning, + "LP solver residuals: primal = %g; dual = %g yield num/max/sum primal " + "(%d/%g/%g) and dual (%d/%g/%g) corrections\n", + norm_primal_residual, norm_dual_residual, int(num_primal_correction), + max_primal_correction, sum_primal_correction, int(num_dual_correction), + max_dual_correction, sum_dual_correction); +} + void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object) { solver_object.model_status_ = HighsModelStatus::kNotset; solver_object.highs_info_.objective_function_value = 0; diff --git a/src/lp_data/HighsSolution.h b/src/lp_data/HighsSolution.h index f869b6b0ee..36342274eb 100644 --- a/src/lp_data/HighsSolution.h +++ b/src/lp_data/HighsSolution.h @@ -123,6 +123,8 @@ HighsStatus formSimplexLpBasisAndFactor( void accommodateAlienBasis(HighsLpSolverObject& solver_object); +void correctResiduals(HighsLpSolverObject& solver_object); + void resetModelStatusAndHighsInfo(HighsLpSolverObject& solver_object); void resetModelStatusAndHighsInfo(HighsModelStatus& model_status, HighsInfo& highs_info); diff --git a/src/lp_data/HighsSolve.cpp b/src/lp_data/HighsSolve.cpp index c03837ac01..9c013b4f2e 100644 --- a/src/lp_data/HighsSolve.cpp +++ b/src/lp_data/HighsSolve.cpp @@ -72,6 +72,11 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { return_status, "solveLpCupdlp"); } if (return_status == HighsStatus::kError) return return_status; + if (solver_object.model_status_ == HighsModelStatus::kOptimal) + // IPM (and PDLP?) can claim optimality with large primal and/or + // dual residual errors, so correct any residual errors that + // exceed the tolerance + correctResiduals(solver_object); // Non-error return requires a primal solution assert(solver_object.solution_.value_valid); // Get the objective and any KKT failures @@ -79,6 +84,10 @@ HighsStatus solveLp(HighsLpSolverObject& solver_object, const string message) { solver_object.lp_.objectiveValue(solver_object.solution_.col_value); getLpKktFailures(options, solver_object.lp_, solver_object.solution_, solver_object.basis_, solver_object.highs_info_); + if (solver_object.model_status_ == HighsModelStatus::kOptimal && + (solver_object.highs_info_.num_primal_infeasibilities > 0 || + solver_object.highs_info_.num_dual_infeasibilities)) + solver_object.model_status_ = HighsModelStatus::kUnknown; if (options.solver == kIpmString || options.run_centring) { // Setting the IPM-specific values of (highs_)info_ has been done in // solveLpIpx diff --git a/src/pdlp/CupdlpWrapper.cpp b/src/pdlp/CupdlpWrapper.cpp index d41f14e595..31fa31e1e8 100644 --- a/src/pdlp/CupdlpWrapper.cpp +++ b/src/pdlp/CupdlpWrapper.cpp @@ -207,7 +207,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, #if CUPDLP_DEBUG analysePdlpSolution(options, lp, highs_solution); #endif - + free(cost); free(lower); free(upper); @@ -225,7 +225,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(prob->lower); free(prob->upper); free(prob->rhs); - + free(prob->hasLower); free(prob->hasUpper); @@ -233,7 +233,7 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(prob->data->csr_matrix->rowMatIdx); free(prob->data->csr_matrix->rowMatElem); free(prob->data->csr_matrix); - + free(prob->data->csc_matrix->colMatBeg); free(prob->data->csc_matrix->colMatIdx); free(prob->data->csc_matrix->colMatElem); @@ -246,13 +246,11 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, free(csc_cpu->colMatBeg); free(csc_cpu->colMatIdx); free(csc_cpu->colMatElem); - + free(csc_cpu); - if (scaling->rowScale != nullptr) - free(scaling->rowScale); - if (scaling->colScale != nullptr) - free(scaling->colScale); + if (scaling->rowScale != nullptr) free(scaling->rowScale); + if (scaling->colScale != nullptr) free(scaling->colScale); free(scaling); return HighsStatus::kOk; diff --git a/src/util/HighsSparseMatrix.cpp b/src/util/HighsSparseMatrix.cpp index 16648e1506..a493d2e7ad 100644 --- a/src/util/HighsSparseMatrix.cpp +++ b/src/util/HighsSparseMatrix.cpp @@ -1174,6 +1174,32 @@ void HighsSparseMatrix::productQuad(vector& result, } } +void HighsSparseMatrix::productTransposeQuad( + vector& result, const vector& row, + const HighsInt debug_report) const { + assert(this->formatOk()); + assert((int)row.size() >= this->num_row_); + result.assign(this->num_col_, 0.0); + if (this->isColwise()) { + for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) { + HighsCDouble value = 0.0; + for (HighsInt iEl = this->start_[iCol]; iEl < this->start_[iCol + 1]; + iEl++) + value += row[this->index_[iEl]] * this->value_[iEl]; + result[iCol] = double(value); + } + } else { + std::vector value(this->num_col_, 0); + for (HighsInt iRow = 0; iRow < this->num_row_; iRow++) { + for (HighsInt iEl = this->start_[iRow]; iEl < this->start_[iRow + 1]; + iEl++) + value[this->index_[iEl]] += row[iRow] * this->value_[iEl]; + } + for (HighsInt iCol = 0; iCol < this->num_col_; iCol++) + result[iCol] = double(value[iCol]); + } +} + void HighsSparseMatrix::productTransposeQuad( vector& result_value, vector& result_index, const HVector& column, const HighsInt debug_report) const { diff --git a/src/util/HighsSparseMatrix.h b/src/util/HighsSparseMatrix.h index 22fef9ea40..7333a37e64 100644 --- a/src/util/HighsSparseMatrix.h +++ b/src/util/HighsSparseMatrix.h @@ -94,6 +94,9 @@ class HighsSparseMatrix { void productTranspose(vector& result, const vector& x) const; void productQuad(vector& result, const vector& x, const HighsInt debug_report = kDebugReportOff) const; + void productTransposeQuad( + vector& result_value, const vector& x, + const HighsInt debug_report = kDebugReportOff) const; void productTransposeQuad( vector& result_value, vector& result_index, const HVector& x, const HighsInt debug_report = kDebugReportOff) const;