Ready to do some testing to identify limit on KKT iterations

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",

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;

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);

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 */

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(); }
 

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; }
 };

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();

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_; }
 

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();

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_;

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) {

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);
   }
 

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(

src/mip/HighsMipSolver.cpp

@@ -806,21 +806,21 @@ 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,
                (long long unsigned)mipdata_->total_lp_iterations,
                (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();

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<double>& 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);

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),

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;