diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp index ed50953852..db48e9c7b7 100644 --- a/src/ipm/IpxWrapper.cpp +++ b/src/ipm/IpxWrapper.cpp @@ -184,6 +184,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer, if (report_solve_data) reportSolveData(options.log_options, ipx_info); highs_info.ipm_iteration_count += (HighsInt)ipx_info.iter; highs_info.crossover_iteration_count += (HighsInt)ipx_info.updates_crossover; + highs_info.max_cr_iteration_count1 = ipx_info.kkt_iter_max1; + highs_info.max_cr_iteration_count2 = ipx_info.kkt_iter_max2; // If not solved... if (solve_status != IPX_STATUS_solved) { @@ -961,6 +963,10 @@ void reportSolveData(const HighsLogOptions& log_options, (int)ipx_info.kktiter1); highsLogDev(log_options, HighsLogType::kInfo, " KKT iter 2 = %d\n", (int)ipx_info.kktiter2); + highsLogDev(log_options, HighsLogType::kInfo, " KKT iter max 1 = %d\n", + (int)ipx_info.kkt_iter_max1); + highsLogDev(log_options, HighsLogType::kInfo, " KKT iter max 2 = %d\n", + (int)ipx_info.kkt_iter_max2); highsLogDev(log_options, HighsLogType::kInfo, " Basis repairs = %d\n", (int)ipx_info.basis_repairs); highsLogDev(log_options, HighsLogType::kInfo, " Updates start = %d\n", diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc index 907f1a156b..cb16878743 100644 --- a/src/ipm/ipx/conjugate_residuals.cc +++ b/src/ipm/ipx/conjugate_residuals.cc @@ -26,13 +26,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, if (maxiter < 0) maxiter = m+100; - Int default_maxiter = maxiter; - Int large_iter = std::max(1000, maxiter/100); - maxiter = std::min(10*large_iter, m+100); - Int report_frequency = 100; - Int report_iter = 0; - printf("ConjugateResiduals::Solve Default maxiter = %d; using %d\n", int(default_maxiter), int(maxiter)); - // Initialize residual, step and Cstep. if (Infnorm(lhs) == 0.0) { residual = rhs; // saves a matrix-vector op @@ -66,12 +59,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs, errflag_ = IPX_ERROR_cr_matrix_not_posdef; break; } - if (iter_ >= report_iter + report_frequency) { - printf("ConjugateResiduals::Solve Iteration count (large = %d; max = %d) is %d: time = %g\n", int(large_iter), int(maxiter), int(iter_), control_.Elapsed()); - report_iter = iter_; - report_frequency *= 2; - } - // Update lhs, residual and Cresidual. const double denom = Dot(Cstep,Cstep); const double alpha = cdot/denom; diff --git a/src/ipm/ipx/info.cc b/src/ipm/ipx/info.cc index fa09ab6cff..a64dea6a58 100644 --- a/src/ipm/ipx/info.cc +++ b/src/ipm/ipx/info.cc @@ -53,6 +53,8 @@ std::ostream& operator<<(std::ostream& os, const Info& info) { dump(os, "iter", info.iter); dump(os, "kktiter1", info.kktiter1); dump(os, "kktiter2", info.kktiter2); + dump(os, "kkt_iter_max1", info.kkt_iter_max1); + dump(os, "kkt_iter_max2", info.kkt_iter_max2); dump(os, "basis_repairs", info.basis_repairs); dump(os, "updates_start", info.updates_start); dump(os, "updates_ipm", info.updates_ipm); diff --git a/src/ipm/ipx/ipx_info.h b/src/ipm/ipx/ipx_info.h index ed43d6718e..03cda96090 100644 --- a/src/ipm/ipx/ipx_info.h +++ b/src/ipm/ipx/ipx_info.h @@ -58,6 +58,8 @@ struct ipx_info { ipxint iter; /* # interior point iterations */ ipxint kktiter1; /* # linear solver iterations before switch */ ipxint kktiter2; /* # linear solver iterations after switch */ + ipxint kkt_iter_max1; /* # max linear solver iterations before switch */ + ipxint kkt_iter_max2; /* # max linear solver iterations after switch */ ipxint basis_repairs; /* # basis repairs after crash, < 0 discarded */ ipxint updates_start; /* # basis updates for starting basis */ ipxint updates_ipm; /* # basis updates in IPM */ diff --git a/src/ipm/ipx/kkt_solver.cc b/src/ipm/ipx/kkt_solver.cc index 1aaee27133..a828fcb0a3 100644 --- a/src/ipm/ipx/kkt_solver.cc +++ b/src/ipm/ipx/kkt_solver.cc @@ -17,7 +17,7 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol, } Int KKTSolver::iterSum() const { return _iterSum(); } - //Int KKTSolver::iterMax() const { return _iterMax(); } +Int KKTSolver::iterMax() const { return _iterMax(); } Int KKTSolver::basis_changes() const { return _basis_changes(); } const Basis* KKTSolver::basis() const { return _basis(); } diff --git a/src/ipm/ipx/kkt_solver.h b/src/ipm/ipx/kkt_solver.h index 5de82440db..e5ab15bc22 100644 --- a/src/ipm/ipx/kkt_solver.h +++ b/src/ipm/ipx/kkt_solver.h @@ -49,9 +49,8 @@ class KKTSolver { Int iterSum() const; // If an iterative method is used, returns the max # iterations in - // all Solve() calls since the last call to Factorize(). A direct - // solver returns the max # iterative refinement steps. - // Int iterMax() const; + // _all_ Solve() calls. + Int iterMax() const; // If a basis matrix is maintained, returns the # basis changes in the last // call to Factorize(). Otherwise returns 0. @@ -66,6 +65,7 @@ class KKTSolver { virtual void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) = 0; virtual Int _iterSum() const = 0; + virtual Int _iterMax() const = 0; virtual Int _basis_changes() const { return 0; } virtual const Basis* _basis() const { return nullptr; } }; diff --git a/src/ipm/ipx/kkt_solver_basis.cc b/src/ipm/ipx/kkt_solver_basis.cc index e335ec9dba..46d28e2766 100644 --- a/src/ipm/ipx/kkt_solver_basis.cc +++ b/src/ipm/ipx/kkt_solver_basis.cc @@ -21,7 +21,6 @@ void KKTSolverBasis::_Factorize(Iterate* iterate, Info* info) { info->errflag = 0; factorized_ = false; iter_sum_ = 0; - iter_max_ = 0; basis_changes_ = 0; for (Int j = 0; j < n+m; j++) @@ -149,6 +148,7 @@ void KKTSolverBasis::_Solve(const Vector& a, const Vector& b, double tol, cr.Solve(splitted_normal_matrix_, work, tol, nullptr, maxiter_, lhs); info->errflag = cr.errflag(); info->kktiter2 += cr.iter(); + info->kkt_iter_max2 = std::max(cr.iter(), info->kkt_iter_max2); info->time_cr2 += cr.time(); info->time_cr2_NNt += splitted_normal_matrix_.time_NNt(); info->time_cr2_B += splitted_normal_matrix_.time_B(); diff --git a/src/ipm/ipx/kkt_solver_basis.h b/src/ipm/ipx/kkt_solver_basis.h index da637324bd..b46a89e5c2 100644 --- a/src/ipm/ipx/kkt_solver_basis.h +++ b/src/ipm/ipx/kkt_solver_basis.h @@ -34,7 +34,7 @@ class KKTSolverBasis : public KKTSolver { void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; Int _iterSum() const override { return iter_sum_; } - // Int _iterMax() const override { return iter_max_; } + Int _iterMax() const override { return iter_max_; } Int _basis_changes() const override { return basis_changes_; } const Basis* _basis() const override { return &basis_; } diff --git a/src/ipm/ipx/kkt_solver_diag.cc b/src/ipm/ipx/kkt_solver_diag.cc index 4ff530cb37..9e66ad6848 100644 --- a/src/ipm/ipx/kkt_solver_diag.cc +++ b/src/ipm/ipx/kkt_solver_diag.cc @@ -19,7 +19,6 @@ void KKTSolverDiag::_Factorize(Iterate* pt, Info* info) { const Int m = model_.rows(); const Int n = model_.cols(); iter_sum_ = 0; - iter_max_ = 0; factorized_ = false; if (pt) { @@ -100,6 +99,7 @@ void KKTSolverDiag::_Solve(const Vector& a, const Vector& b, double tol, cr.Solve(normal_matrix_, precond_, rhs, tol, &resscale_[0], maxiter_, y); info->errflag = cr.errflag(); info->kktiter1 += cr.iter(); + info->kkt_iter_max1 = std::max(cr.iter(), info->kkt_iter_max1); info->time_cr1 += cr.time(); info->time_cr1_AAt += normal_matrix_.time(); info->time_cr1_pre += precond_.time(); diff --git a/src/ipm/ipx/kkt_solver_diag.h b/src/ipm/ipx/kkt_solver_diag.h index 66029e0a6a..aaa18cfb85 100644 --- a/src/ipm/ipx/kkt_solver_diag.h +++ b/src/ipm/ipx/kkt_solver_diag.h @@ -30,7 +30,7 @@ class KKTSolverDiag : public KKTSolver { void _Solve(const Vector& a, const Vector& b, double tol, Vector& x, Vector& y, Info* info) override; Int _iterSum() const override { return iter_sum_; }; - // Int _iterMax() const override { return iter_max_; }; + Int _iterMax() const override { return iter_max_; }; const Control& control_; const Model& model_; diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp index 987c18cc70..bbb07eb170 100644 --- a/src/lp_data/HighsInfo.cpp +++ b/src/lp_data/HighsInfo.cpp @@ -41,7 +41,8 @@ void HighsInfo::invalidate() { max_complementarity_violation = kHighsIllegalComplementarityViolation; sum_complementarity_violations = kHighsIllegalComplementarityViolation; primal_dual_integral = -kHighsInf; - max_cr_iteration_count = -1; + max_cr_iteration_count1 = -1; + max_cr_iteration_count2 = -1; } static std::string infoEntryTypeToString(const HighsInfoType type) { diff --git a/src/lp_data/HighsInfo.h b/src/lp_data/HighsInfo.h index 294b4350ed..70b56c20ef 100644 --- a/src/lp_data/HighsInfo.h +++ b/src/lp_data/HighsInfo.h @@ -116,7 +116,8 @@ struct HighsInfoStruct { double max_complementarity_violation; double sum_complementarity_violations; double primal_dual_integral; - HighsInt max_cr_iteration_count; + HighsInt max_cr_iteration_count1; + HighsInt max_cr_iteration_count2; }; class HighsInfo : public HighsInfoStruct { @@ -281,9 +282,16 @@ class HighsInfo : public HighsInfoStruct { advanced, &primal_dual_integral, 0); records.push_back(record_double); - record_int = new InfoRecordInt("max_cr_iteration_count", - "Maximum number of CG iterations in IPX", advanced, - &max_cr_iteration_count, -1); + record_int = + new InfoRecordInt("max_cr_iteration_count1", + "Maximum number of diag CR iterations in IPX", + advanced, &max_cr_iteration_count1, -1); + records.push_back(record_int); + + record_int = + new InfoRecordInt("max_cr_iteration_count2", + "Maximum number of basic CR iterations in IPX", + advanced, &max_cr_iteration_count2, -1); records.push_back(record_int); } diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h index 7976c0bd8c..73c261a889 100644 --- a/src/lp_data/HighsOptions.h +++ b/src/lp_data/HighsOptions.h @@ -558,7 +558,7 @@ struct HighsOptionsStruct { mip_pool_soft_limit(0), mip_pscost_minreliable(0), mip_min_cliquetable_entries_for_parallelism(0), - mip_compute_analytic_centre(0), + mip_compute_analytic_centre(0), mip_report_level(0), mip_feasibility_tolerance(0.0), mip_rel_gap(0.0), @@ -1038,11 +1038,10 @@ class HighsOptions : public HighsOptionsStruct { record_int = new OptionRecordInt( "mip_compute_analytic_centre", - "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 => true", - advanced, &mip_compute_analytic_centre, - kMipAnalyticCentreCalulationMin, - kMipAnalyticCentreCalulationOriginal, - kMipAnalyticCentreCalulationMax); + "Compute analytic centre for MIP: 0 => no; 1 => original (default) 2 " + "=> true", + advanced, &mip_compute_analytic_centre, kMipAnalyticCentreCalulationMin, + kMipAnalyticCentreCalulationOriginal, kMipAnalyticCentreCalulationMax); records.push_back(record_int); record_int = @@ -1084,9 +1083,11 @@ class HighsOptions : public HighsOptionsStruct { &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf); records.push_back(record_int); - record_int = new OptionRecordInt( - "kkt_iteration_limit", "Iteration limit for PCG in IPX IPM solver: default is -1, so set internally", advanced, - &kkt_iteration_limit, -1, -1, kHighsIInf); + record_int = + new OptionRecordInt("kkt_iteration_limit", + "Iteration limit for PCG in IPX IPM solver: " + "default is -1, so set internally", + advanced, &kkt_iteration_limit, -1, -1, kHighsIInf); records.push_back(record_int); record_bool = new OptionRecordBool( diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp index ead83e320a..66520a4470 100644 --- a/src/mip/HighsMipSolver.cpp +++ b/src/mip/HighsMipSolver.cpp @@ -806,13 +806,13 @@ void HighsMipSolver::cleanupSolve() { " %llu (strong br.)\n" " %llu (separation)\n" " %llu (heuristics)\n" - " AC failures %d\n" + " AC failures %d\n" " Max sub-MIP depth %d\n", timer_.read(timer_.total_clock), analysis_.mipTimerRead(kMipClockPresolve), analysis_.mipTimerRead(kMipClockSolve), analysis_.mipTimerRead(kMipClockPostsolve), - (long long unsigned)mipdata_->num_nodes, + (long long unsigned)mipdata_->num_nodes, (long long unsigned)mipdata_->total_repair_lp, (long long unsigned)mipdata_->total_repair_lp_feasible, (long long unsigned)mipdata_->total_repair_lp_iterations, @@ -820,7 +820,7 @@ 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(mipdata_->num_analytic_centre_failures), int(max_submip_level > 0 ? max_submip_level : 0)); analysis_.reportMipTimer(); diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp index f7438d7358..5fe9bff4b9 100644 --- a/src/mip/HighsMipSolverData.cpp +++ b/src/mip/HighsMipSolverData.cpp @@ -302,9 +302,11 @@ void HighsMipSolverData::startAnalyticCenterComputation( taskGroup.spawn([&]() { // first check if the analytic centre computation should be cancelled, e.g. // due to early return in the root node evaluation - assert(mipsolver.options_mip_->mip_compute_analytic_centre >= kMipAnalyticCentreCalulationOriginal); + assert(mipsolver.options_mip_->mip_compute_analytic_centre >= + kMipAnalyticCentreCalulationOriginal); Highs ipm; - if (mipsolver.options_mip_->mip_compute_analytic_centre == kMipAnalyticCentreCalulationOriginal) { + if (mipsolver.options_mip_->mip_compute_analytic_centre == + kMipAnalyticCentreCalulationOriginal) { // Original calculation is just IPM with crossover off ipm.setOptionValue("solver", "ipm"); ipm.setOptionValue("run_crossover", kHighsOffString); @@ -321,11 +323,18 @@ void HighsMipSolverData::startAnalyticCenterComputation( lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0); ipm.passModel(std::move(lpmodel)); - // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + "_ipm.mps"); + // if (!mipsolver.submip) ipm.writeModel(mipsolver.model_->model_name_ + + // "_ipm.mps"); 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)); + } const std::vector& sol = ipm.getSolution().col_value; if (HighsInt(sol.size()) != mipsolver.numCol()) return; analyticCenterStatus = ipm.getModelStatus(); @@ -1877,9 +1886,10 @@ void HighsMipSolverData::evaluateRootNode() { if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed && !analyticCenterFailed) { if (mipsolver.analysis_.analyse_mip_time) { - highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo, - "MIP-Timing: %11.2g - starting analytic centre calculation\n", - mipsolver.timer_.read()); + highsLogUser( + mipsolver.options_mip_->log_options, HighsLogType::kInfo, + "MIP-Timing: %11.2g - starting analytic centre calculation\n", + mipsolver.timer_.read()); fflush(stdout); } mipsolver.analysis_.mipTimerStart(kMipClockStartAnalyticCentreComputation); @@ -2059,7 +2069,8 @@ void HighsMipSolverData::evaluateRootNode() { return; } if (nseparounds >= 5 && !mipsolver.submip && - mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { if (checkLimits()) { mipsolver.analysis_.mipTimerStop(kMipClockSeparation); return; @@ -2154,7 +2165,8 @@ void HighsMipSolverData::evaluateRootNode() { rootlpsolobj = lp.getObjective(); lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters))); - if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { if (checkLimits()) return; mipsolver.analysis_.mipTimerStart(kMipClockFinishAnalyticCentreComputation); @@ -2276,7 +2288,8 @@ void HighsMipSolverData::evaluateRootNode() { if (lower_bound <= upper_limit) { if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart && mipsolver.options_mip_->presolve != kHighsOffString) { - if (mipsolver.options_mip_->mip_compute_analytic_centre && !analyticCenterComputed) { + if (mipsolver.options_mip_->mip_compute_analytic_centre && + !analyticCenterComputed) { mipsolver.analysis_.mipTimerStart( kMipClockFinishAnalyticCentreComputation); finishAnalyticCenterComputation(tg); diff --git a/src/mip/HighsMipSolverData.h b/src/mip/HighsMipSolverData.h index 6e101af638..915d7c9761 100644 --- a/src/mip/HighsMipSolverData.h +++ b/src/mip/HighsMipSolverData.h @@ -172,7 +172,7 @@ struct HighsMipSolverData { cliquesExtracted(false), rowMatrixSet(false), analyticCenterComputed(false), - analyticCenterFailed(false), + analyticCenterFailed(false), analyticCenterStatus(HighsModelStatus::kNotset), detectSymmetries(false), numRestarts(0), @@ -187,7 +187,7 @@ struct HighsMipSolverData { rootlpsolobj(-kHighsInf), numintegercols(0), maxTreeSizeLog2(0), - num_analytic_centre_failures(0), + num_analytic_centre_failures(0), pruned_treeweight(0), avgrootlpiters(0.0), last_disptime(0.0), diff --git a/src/mip/HighsPrimalHeuristics.cpp b/src/mip/HighsPrimalHeuristics.cpp index 3729858ec1..d3a90ef8b7 100644 --- a/src/mip/HighsPrimalHeuristics.cpp +++ b/src/mip/HighsPrimalHeuristics.cpp @@ -135,7 +135,8 @@ bool HighsPrimalHeuristics::solveSubMip( // If the analytic centre calculation has failed for the parent MIP, // then don't perform it for the sub-MIP - // if (mipsolver.mipdata_->analyticCenterFailed) submipoptions.mip_compute_analytic_centre = 0; + // if (mipsolver.mipdata_->analyticCenterFailed) + // submipoptions.mip_compute_analytic_centre = 0; // setup solver and run it HighsSolution solution; solution.value_valid = false;