diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index cb16878743..f7345def7c 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -9,6 +9,8 @@ namespace ipx { ConjugateResiduals::ConjugateResiduals(const Control& control) : control_(control) {} +const bool cr_logging = false; + void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, double tol, const double* resscale, Int maxiter, Vector& lhs) { @@ -23,9 +25,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = 0; iter_ = 0; time_ = 0.0; + Int use_maxiter = maxiter; if (maxiter < 0) - maxiter = m+100; + use_maxiter = m+100; + if (cr_logging) printf("\nCR(C) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout); + maxiter = use_maxiter; // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -37,9 +42,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, step = residual; Cstep = Cresidual; + double resnorm; while (true) { // Termination check. - double resnorm = 0.0; + resnorm = 0.0; if (resscale) for (Int i = 0; i < m; i++) resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i])); @@ -48,6 +54,7 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (resnorm <= tol) break; if (iter_ == maxiter) { + if (cr_logging) printf("CR(C) reached maxiter!\n\n"); fflush(stdout); control_.Debug(3) << " CR method not converged in " << maxiter << " iterations." << " residual = " << sci2(resnorm) << ',' @@ -78,11 +85,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, cdot = cdotnew; iter_++; - if ((errflag_ = control_.InterruptCheck()) != 0) { - printf("ConjugateResiduals::Solve control_.InterruptCheck()) != 0: time = %g\n", control_.Elapsed()); - break; - } + if ((errflag_ = control_.InterruptCheck()) != 0) break; } + if (errflag_ != IPX_ERROR_cr_iter_limit) + if (cr_logging) printf("CR(C) iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout); time_ = timer.Elapsed(); } @@ -110,9 +116,12 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, errflag_ = 0; iter_ = 0; time_ = 0.0; + Int use_maxiter = maxiter; if (maxiter < 0) - maxiter = m+100; + use_maxiter = m+100; + if (cr_logging) printf("\nCR(C,P) maxiter = (entry = %d, use = %d)\n", int(maxiter), int(use_maxiter)); fflush(stdout); + maxiter = use_maxiter; // Initialize residual, sresidual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -125,9 +134,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, step = sresidual; Cstep = Csresidual; + double resnorm; while (true) { // Termination check. - double resnorm = 0.0; + resnorm = 0.0; if (resscale) for (Int i = 0; i < m; i++) resnorm = std::max(resnorm, std::abs(resscale[i]*residual[i])); @@ -136,7 +146,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, if (resnorm <= tol) break; if (iter_ == maxiter) { - control_.Debug(3) + if (cr_logging) printf("CR(C,P) reached maxiter!\n\n"); fflush(stdout); + control_.Debug(3) << " PCR method not converged in " << maxiter << " iterations." << " residual = " << sci2(resnorm) << ',' << " tolerance = " << sci2(tol) << '\n'; @@ -208,6 +219,8 @@ void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P, if ((errflag_ = control_.InterruptCheck()) != 0) break; } + if (errflag_ != IPX_ERROR_cr_iter_limit) + if (cr_logging) printf("CR(C,P) iter = %d; resnorm = %g <= %g = tol? %s\n\n", int(iter_), resnorm, tol, resnorm <= tol ? "T" : "F"); fflush(stdout); time_ = timer.Elapsed(); } diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index 66520a4470..6585c0ddb5 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -795,6 +795,11 @@ void HighsMipSolver::cleanupSolve() { " %.12g (row viol.)\n", solution_objective_, bound_violation_, integrality_violation_, row_violation_); + int ac_start = mipdata_->num_analytic_centre_start; + int ac_opt = mipdata_->num_analytic_centre_opt; + int ac_other = mipdata_->num_analytic_centre_other; + int ac_fail = mipdata_->num_analytic_centre_fail; + int ac_abort = ac_start - ac_opt - ac_other - ac_fail; highsLogUser(options_mip_->log_options, HighsLogType::kInfo, " Timing %.2f (total)\n" " %.2f (presolve)\n" @@ -806,8 +811,9 @@ void HighsMipSolver::cleanupSolve() { " %llu (strong br.)\n" " %llu (separation)\n" " %llu (heuristics)\n" - " AC failures %d\n" - " Max sub-MIP depth %d\n", + " Max sub-MIP depth %d\n" + " Analytic centre calculations (Optimal / other / abort / " + "fail) = (%d / %d / %d / %d) of %d\n", timer_.read(timer_.total_clock), analysis_.mipTimerRead(kMipClockPresolve), analysis_.mipTimerRead(kMipClockSolve), @@ -820,8 +826,12 @@ void HighsMipSolver::cleanupSolve() { (long long unsigned)mipdata_->sb_lp_iterations, (long long unsigned)mipdata_->sepa_lp_iterations, (long long unsigned)mipdata_->heuristic_lp_iterations, - int(mipdata_->num_analytic_centre_failures), - int(max_submip_level > 0 ? max_submip_level : 0)); + int(max_submip_level > 0 ? max_submip_level : 0), ac_opt, + ac_other, ac_abort, ac_fail, ac_start); + + highsLogUser(options_mip_->log_options, HighsLogType::kInfo, + "grepAcStats,%s,%d,%d,%d,%d,%d\n", model_->model_name_.c_str(), + ac_opt, ac_other, ac_abort, ac_fail, ac_start); analysis_.reportMipTimer(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index 5fe9bff4b9..f4af08b640 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -304,6 +304,7 @@ void HighsMipSolverData::startAnalyticCenterComputation( // due to early return in the root node evaluation assert(mipsolver.options_mip_->mip_compute_analytic_centre >= kMipAnalyticCentreCalulationOriginal); + const bool ac_logging = !mipsolver.submip; Highs ipm; if (mipsolver.options_mip_->mip_compute_analytic_centre == kMipAnalyticCentreCalulationOriginal) { @@ -317,8 +318,11 @@ void HighsMipSolverData::startAnalyticCenterComputation( } ipm.setOptionValue("presolve", "off"); ipm.setOptionValue("output_flag", false); - ipm.setOptionValue("output_flag", !mipsolver.submip); + // ipm.setOptionValue("output_flag", !mipsolver.submip); ipm.setOptionValue("ipm_iteration_limit", 200); + HighsInt kkt_iteration_limit = + std::max(100, mipsolver.model_->num_row_ / 1000); + ipm.setOptionValue("kkt_iteration_limit", kkt_iteration_limit); HighsLp lpmodel(*mipsolver.model_); lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); @@ -326,14 +330,28 @@ void HighsMipSolverData::startAnalyticCenterComputation( // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + // "_ipm.mps"); + if (ac_logging) num_analytic_centre_start++; + double tt0 = mipsolver.timer_.read(mipsolver.timer_.total_clock); mipsolver.analysis_.mipTimerStart(kMipClockIpmSolveLp); ipm.run(); mipsolver.analysis_.mipTimerStop(kMipClockIpmSolveLp); - if (!mipsolver.submip) { - printf("grepAC: model; num_row; max CR counts are, %s, %d, %d, %d\n", - mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), - int(ipm.getInfo().max_cr_iteration_count1), - int(ipm.getInfo().max_cr_iteration_count2)); + double tt1 = mipsolver.timer_.read(mipsolver.timer_.total_clock); + if (ac_logging) { + HighsModelStatus ac_status = ipm.getModelStatus(); + printf( + "grepAcLocal: model; num_row; CR limit; max CR counts; time and " + "status, %s, %d, %d, %d, %d, %g, %s\n", + mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()), + int(kkt_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1), + int(ipm.getInfo().max_cr_iteration_count2), tt1 - tt0, + ipm.modelStatusToString(ac_status).c_str()); + if (ac_status == HighsModelStatus::kSolveError) { + num_analytic_centre_fail++; + } else if (ac_status == HighsModelStatus::kOptimal) { + num_analytic_centre_opt++; + } else { + num_analytic_centre_other++; + } } const std::vector& sol = ipm.getSolution().col_value; if (HighsInt(sol.size()) != mipsolver.numCol()) return; @@ -359,7 +377,6 @@ void HighsMipSolverData::finishAnalyticCenterComputation( } analyticCenterComputed = true; analyticCenterFailed = analyticCenterStatus == HighsModelStatus::kSolveError; - if (analyticCenterFailed) num_analytic_centre_failures++; if (analyticCenterStatus == HighsModelStatus::kOptimal) { HighsInt nfixed = 0; HighsInt nintfixed = 0; @@ -649,7 +666,10 @@ void HighsMipSolverData::init() { heuristic_lp_iterations_before_run = 0; sepa_lp_iterations_before_run = 0; sb_lp_iterations_before_run = 0; - num_analytic_centre_failures = 0; + num_analytic_centre_start = 0; + num_analytic_centre_opt = 0; + num_analytic_centre_other = 0; + num_analytic_centre_fail = 0; num_disp_lines = 0; numCliqueEntriesAfterPresolve = 0; numCliqueEntriesAfterFirstPresolve = 0; diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 915d7c9761..ba8b8c3323 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -120,7 +120,10 @@ struct HighsMipSolverData { HighsInt numintegercols; HighsInt maxTreeSizeLog2; - HighsInt num_analytic_centre_failures; + HighsInt num_analytic_centre_start; + HighsInt num_analytic_centre_opt; + HighsInt num_analytic_centre_other; + HighsInt num_analytic_centre_fail; HighsCDouble pruned_treeweight; double avgrootlpiters; double last_disptime; @@ -187,7 +190,10 @@ struct HighsMipSolverData { rootlpsolobj(-kHighsInf), numintegercols(0), maxTreeSizeLog2(0), - num_analytic_centre_failures(0), + num_analytic_centre_start(0), + num_analytic_centre_opt(0), + num_analytic_centre_other(0), + num_analytic_centre_fail(0), pruned_treeweight(0), avgrootlpiters(0.0), last_disptime(0.0), diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index d3a90ef8b7..91f53a0ea2 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -137,6 +137,7 @@ bool HighsPrimalHeuristics::solveSubMip( // if (mipsolver.mipdata_->analyticCenterFailed) // submipoptions.mip_compute_analytic_centre = 0; + // setup solver and run it HighsSolution solution; solution.value_valid = false;