Introducing a deterministic decision on whether to compute the analytic centre in MIP

src/ipm/IpxWrapper.cpp

@@ -121,6 +121,7 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   parameters.time_limit = options.time_limit - timer.readRunHighsClock();
   parameters.ipm_maxiter =
       options.ipm_iteration_limit - highs_info.ipm_iteration_count;
+  parameters.kkt_maxiter = options.kkt_iteration_limit;
   // Determine if crossover is to be run or not
   //
   // When doing analytic centring calculations, crossover must not be
@@ -183,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) {
@@ -960,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

@@ -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) << ','
@@ -59,7 +66,6 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
             errflag_ = IPX_ERROR_cr_matrix_not_posdef;
             break;
         }
-
         // Update lhs, residual and Cresidual.
         const double denom = Dot(Cstep,Cstep);
         const double alpha = cdot/denom;
@@ -79,9 +85,10 @@ void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
         cdot = cdotnew;
 
         iter_++;
-        if ((errflag_ = control_.InterruptCheck()) != 0)
-            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();
 }
 
@@ -109,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
@@ -124,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]));
@@ -135,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';
@@ -207,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();
 }
 

src/ipm/ipx/control.cc

@@ -13,8 +13,10 @@ Control::Control() {
 Int Control::InterruptCheck(const Int ipm_iteration_count) const {
     HighsTaskExecutor::getThisWorkerDeque()->checkInterrupt();
     if (parameters_.time_limit >= 0.0 &&
-        parameters_.time_limit < timer_.Elapsed())
+        parameters_.time_limit < timer_.Elapsed()) {
+      printf("Control::InterruptCheck Reached time limit of %g\n", parameters_.time_limit);
         return IPX_ERROR_time_interrupt;
+    }
     // The pointer callback_ should not be null, since that indicates
     // that it's not been set
     assert(callback_);

src/ipm/ipx/control.h

@@ -69,6 +69,7 @@ class Control {
     double ipm_drop_primal() const { return parameters_.ipm_drop_primal; }
     double ipm_drop_dual() const { return parameters_.ipm_drop_dual; }
     double kkt_tol() const { return parameters_.kkt_tol; }
+    ipxint kkt_maxiter() const { return parameters_.kkt_maxiter; }
     ipxint crash_basis() const { return parameters_.crash_basis; }
     double dependency_tol() const { return parameters_.dependency_tol; }
     double volume_tol() const { return parameters_.volume_tol; }

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/ipm.cc

@@ -856,7 +856,7 @@ void IPM::PrintOutput() {
     control_.Debug()
       << "  " << Fixed(step_primal_, 4, 2) << " " << Fixed(step_dual_, 4, 2)
       << "  " << Format(kkt_->basis_changes(), 7)
-      << " "  << Format(kkt_->iter(), 7);
+      << " "  << Format(kkt_->iterSum(), 7);
     control_.Debug()
       << "  " << Format(info_->dual_dropped, 7)
       << " "  << Format(info_->primal_dropped, 7);

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/ipx_internal.h

@@ -32,6 +32,7 @@ struct Parameters : public ipx_parameters {
     ipm_drop_primal = 1e-9;
     ipm_drop_dual = 1e-9;
     kkt_tol = 0.3;
+    kkt_maxiter = -1;
     crash_basis = 1;
     dependency_tol = 1e-6;
     volume_tol = 2.0;

src/ipm/ipx/ipx_parameters.h

@@ -28,6 +28,7 @@ struct ipx_parameters {
     double ipm_drop_dual;
 
     /* Linear solver */
+    ipxint kkt_maxiter;
     double kkt_tol;
 
     /* Basis construction in IPM */

src/ipm/ipx/kkt_solver.cc

@@ -16,7 +16,8 @@ void KKTSolver::Solve(const Vector& a, const Vector& b, double tol,
     info->time_kkt_solve += timer.Elapsed();
 }
 
-Int KKTSolver::iter() const { return _iter(); }
+Int KKTSolver::iterSum() const { return _iterSum(); }
+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

@@ -46,7 +46,11 @@ class KKTSolver {
     // If an iterative method is used, returns the # iterations in all Solve()
     // calls since the last call to Factorize(). A direct solver returns the #
     // iterative refinement steps.
-    Int iter() const;
+    Int iterSum() const;
+
+    // If an iterative method is used, returns the max # iterations in
+    // _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.
@@ -60,7 +64,8 @@ class KKTSolver {
     virtual void _Factorize(Iterate* iterate, Info* info) = 0;
     virtual void _Solve(const Vector& a, const Vector& b, double tol,
                          Vector& x, Vector& y, Info* info) = 0;
-    virtual Int _iter() const = 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

@@ -20,7 +20,7 @@ void KKTSolverBasis::_Factorize(Iterate* iterate, Info* info) {
     const Int n = model_.cols();
     info->errflag = 0;
     factorized_ = false;
-    iter_ = 0;
+    iter_sum_ = 0;
     basis_changes_ = 0;
 
     for (Int j = 0; j < n+m; j++)
@@ -148,11 +148,13 @@ 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();
     info->time_cr2_Bt += splitted_normal_matrix_.time_Bt();
-    iter_ += cr.iter();
+    iter_sum_ += cr.iter();
+    iter_max_ = std::max(cr.iter(), iter_max_);
 
     // Permute back solution to normal equations.
     for (Int k = 0; k < m; k++)

src/ipm/ipx/kkt_solver_basis.h

@@ -33,7 +33,8 @@ class KKTSolverBasis : public KKTSolver {
     void _Factorize(Iterate* iterate, Info* info) override;
     void _Solve(const Vector& a, const Vector& b, double tol,
                 Vector& x, Vector& y, Info* info) override;
-    Int _iter() const override { return iter_; }
+    Int _iterSum() const override { return iter_sum_; }
+    Int _iterMax() const override { return iter_max_; }
     Int _basis_changes() const override { return basis_changes_; }
     const Basis* _basis() const override { return &basis_; }
 
@@ -57,7 +58,8 @@ class KKTSolverBasis : public KKTSolver {
     Vector colscale_;           // interior point column scaling factors
     bool factorized_{false};    // preconditioner prepared?
     Int maxiter_{-1};
-    Int iter_{0};
+    Int iter_sum_{0};
+    Int iter_max_{0};
     Int basis_changes_{0};
 };
 

src/ipm/ipx/kkt_solver_diag.cc

@@ -18,7 +18,7 @@ KKTSolverDiag::KKTSolverDiag(const Control& control, const Model& model) :
 void KKTSolverDiag::_Factorize(Iterate* pt, Info* info) {
     const Int m = model_.rows();
     const Int n = model_.cols();
-    iter_ = 0;
+    iter_sum_ = 0;
     factorized_ = false;
 
     if (pt) {
@@ -99,10 +99,12 @@ 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();
-    iter_ += cr.iter();
+    iter_sum_ += cr.iter();
+    iter_max_ = std::max(cr.iter(), iter_max_);
 
     // Recover solution to KKT system.
     for (Int i = 0; i < m; i++)

src/ipm/ipx/kkt_solver_diag.h

@@ -29,7 +29,8 @@ class KKTSolverDiag : public KKTSolver {
     void _Factorize(Iterate* iterate, Info* info) override;
     void _Solve(const Vector& a, const Vector& b, double tol,
                 Vector& x, Vector& y, Info* info) override;
-    Int _iter() const override { return iter_; };
+    Int _iterSum() const override { return iter_sum_; };
+    Int _iterMax() const override { return iter_max_; };
 
     const Control& control_;
     const Model& model_;
@@ -40,7 +41,8 @@ class KKTSolverDiag : public KKTSolver {
     Vector resscale_;        // residual scaling factors for CR termination test
     bool factorized_{false}; // KKT matrix factorized?
     Int maxiter_{-1};
-    Int iter_{0};               // # CR iterations since last Factorize()
+    Int iter_sum_{0};
+    Int iter_max_{0};
 };
 
 }  // namespace ipx

src/ipm/ipx/lp_solver.cc

@@ -65,16 +65,16 @@ Int LpSolver::Solve() {
 	//        if ((info_.status_ipm == IPX_STATUS_optimal ||
 	//             info_.status_ipm == IPX_STATUS_imprecise) && run_crossover_on) {
 	if (run_crossover) {
-	    if (run_crossover_on) {
-	      control_.hLog("Running crossover as requested\n");
-	    } else if (run_crossover_choose) {
-	      assert(info_.status_ipm == IPX_STATUS_imprecise);
-	      control_.hLog("Running crossover since IPX is imprecise\n");
-	    } else {
-	      assert(run_crossover_on || run_crossover_choose);
-	    }
-	    BuildCrossoverStartingPoint();
-            RunCrossover();
+	  if (run_crossover_on) {
+	    control_.hLog("Running crossover as requested\n");
+	  } else if (run_crossover_choose) {
+	    assert(info_.status_ipm == IPX_STATUS_imprecise);
+	    control_.hLog("Running crossover since IPX is imprecise\n");
+	  } else {
+	    assert(run_crossover_on || run_crossover_choose);
+	  }
+	  BuildCrossoverStartingPoint();
+	  RunCrossover();
         }
         if (basis_) {
             info_.ftran_sparse = basis_->frac_ftran_sparse();
@@ -558,6 +558,7 @@ void LpSolver::RunMainIPM(IPM& ipm) {
     KKTSolverBasis kkt(control_, *basis_);
     Timer timer;
     ipm.maxiter(control_.ipm_maxiter());
+    kkt.maxiter(control_.kkt_maxiter());
     ipm.Driver(&kkt, iterate_.get(), &info_);
     info_.time_ipm2 = timer.Elapsed();
 }

src/lp_data/HConst.h

@@ -45,6 +45,14 @@ const double kExcessivelySmallCostValue = 1e-4;
 const bool kAllowDeveloperAssert = false;
 const bool kExtendInvertWhenAddingRows = false;
 
+enum MipAnalyticCentreCalulation {
+  kMipAnalyticCentreCalulationMin = 0,
+  kMipAnalyticCentreCalulationNo = kMipAnalyticCentreCalulationMin,
+  kMipAnalyticCentreCalulationOriginal,
+  kMipAnalyticCentreCalulationTrue,
+  kMipAnalyticCentreCalulationMax = kMipAnalyticCentreCalulationTrue
+};
+
 enum class HighsLogType { kInfo = 1, kDetailed, kVerbose, kWarning, kError };
 
 enum SimplexScaleStrategy {

src/lp_data/HighsInfo.cpp

@@ -41,6 +41,8 @@ void HighsInfo::invalidate() {
   max_complementarity_violation = kHighsIllegalComplementarityViolation;
   sum_complementarity_violations = kHighsIllegalComplementarityViolation;
   primal_dual_integral = -kHighsInf;
+  max_cr_iteration_count1 = -1;
+  max_cr_iteration_count2 = -1;
 }
 
 static std::string infoEntryTypeToString(const HighsInfoType type) {