diff --git a/check/TestLpSolvers.cpp b/check/TestLpSolvers.cpp
index 49da8322b0..2a59938efe 100644
--- a/check/TestLpSolvers.cpp
+++ b/check/TestLpSolvers.cpp
@@ -2,7 +2,7 @@
 #include "Highs.h"
 #include "catch.hpp"
 
-const bool dev_run = false;
+const bool dev_run = true;
 
 struct IterationCount {
   HighsInt simplex;
@@ -648,27 +648,58 @@ TEST_CASE("simplex-stats", "[highs_lp_solver]") {
 
   Highs h;
   const HighsSimplexStats& simplex_stats = h.getSimplexStats();
+  const HighsSimplexStats& presolved_lp_simplex_stats =
+      h.getPresolvedLpSimplexStats();
   h.setOptionValue("output_flag", dev_run);
+  // std::string model = "dcp2";
+  std::string model = "adlittle";
   std::string model_file =
-      std::string(HIGHS_DIR) + "/check/instances/adlittle.mps";
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
   REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
 
   REQUIRE(h.run() == HighsStatus::kOk);
   REQUIRE(simplex_stats.valid);
-  REQUIRE(simplex_stats.iteration_count == 0);
-  REQUIRE(simplex_stats.num_invert == 1);
+  REQUIRE(simplex_stats.num_col > 0);
+  REQUIRE(simplex_stats.num_row > 0);
+  REQUIRE(simplex_stats.num_nz > 0);
+  REQUIRE(simplex_stats.iteration_count >= 0);
+  REQUIRE(simplex_stats.num_invert > 0);
   REQUIRE(simplex_stats.last_invert_num_el > 0);
   REQUIRE(simplex_stats.last_factored_basis_num_el > 0);
-  REQUIRE(simplex_stats.col_aq_density == 0);
-  REQUIRE(simplex_stats.row_ep_density == 0);
-  REQUIRE(simplex_stats.row_ap_density == 0);
-  REQUIRE(simplex_stats.row_DSE_density == 0);
+  if (simplex_stats.iteration_count > 0) {
+    REQUIRE(simplex_stats.col_aq_density > 0);
+    REQUIRE(simplex_stats.row_ep_density > 0);
+    REQUIRE(simplex_stats.row_ap_density > 0);
+    REQUIRE(simplex_stats.row_DSE_density > 0);
+  } else {
+    REQUIRE(simplex_stats.col_aq_density == 0);
+    REQUIRE(simplex_stats.row_ep_density == 0);
+    REQUIRE(simplex_stats.row_ap_density == 0);
+    REQUIRE(simplex_stats.row_DSE_density == 0);
+  }
   if (dev_run) h.reportSimplexStats(stdout);
 
+  REQUIRE(presolved_lp_simplex_stats.valid);
+  REQUIRE(presolved_lp_simplex_stats.num_col > 0);
+  REQUIRE(presolved_lp_simplex_stats.num_row > 0);
+  REQUIRE(presolved_lp_simplex_stats.num_nz > 0);
+  REQUIRE(presolved_lp_simplex_stats.iteration_count > 0);
+  REQUIRE(presolved_lp_simplex_stats.num_invert > 0);
+  REQUIRE(presolved_lp_simplex_stats.last_invert_num_el > 0);
+  REQUIRE(presolved_lp_simplex_stats.last_factored_basis_num_el > 0);
+  REQUIRE(presolved_lp_simplex_stats.col_aq_density > 0);
+  REQUIRE(presolved_lp_simplex_stats.row_ep_density > 0);
+  REQUIRE(presolved_lp_simplex_stats.row_ap_density > 0);
+  REQUIRE(presolved_lp_simplex_stats.row_DSE_density > 0);
+  if (dev_run) h.reportPresolvedLpSimplexStats(stdout);
+
   h.clearSolver();
   h.setOptionValue("presolve", kHighsOffString);
   REQUIRE(h.run() == HighsStatus::kOk);
   REQUIRE(simplex_stats.valid);
+  REQUIRE(simplex_stats.num_col > 0);
+  REQUIRE(simplex_stats.num_row > 0);
+  REQUIRE(simplex_stats.num_nz > 0);
   REQUIRE(simplex_stats.iteration_count > 0);
   REQUIRE(simplex_stats.num_invert > 0);
   REQUIRE(simplex_stats.last_invert_num_el > 0);
@@ -679,3 +710,36 @@ TEST_CASE("simplex-stats", "[highs_lp_solver]") {
   REQUIRE(simplex_stats.row_DSE_density > 0);
   if (dev_run) h.reportSimplexStats(stdout);
 }
+
+TEST_CASE("ipx-stats", "[highs_lp_solver]") {
+  HighsStatus return_status;
+
+  Highs h;
+  const HighsIpxStats& ipx_stats = h.getIpxStats();
+  h.setOptionValue("output_flag", dev_run);
+  // std::string model = "dcp2";
+  std::string model = "adlittle";
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+  h.setOptionValue("solver", kIpmString);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(ipx_stats.valid);
+  REQUIRE(ipx_stats.num_col > 0);
+  REQUIRE(ipx_stats.num_row > 0);
+  REQUIRE(ipx_stats.num_nz > 0);
+  REQUIRE(ipx_stats.iteration_count > 0);
+  for (HighsInt iteration = 0; iteration < ipx_stats.iteration_count;
+       iteration++) {
+    REQUIRE(ipx_stats.cr_count[iteration] > 0);
+    if (ipx_stats.cr_type[iteration] == 2) {
+      REQUIRE(ipx_stats.factored_basis_num_el[iteration] > 0);
+      REQUIRE(ipx_stats.invert_num_el[iteration] > 0);
+    } else {
+      REQUIRE(ipx_stats.cr_type[iteration] == 1);
+      REQUIRE(ipx_stats.factored_basis_num_el[iteration] == 0);
+      REQUIRE(ipx_stats.invert_num_el[iteration] == 0);
+    }
+  }
+  if (dev_run) h.reportIpxStats(stdout);
+}
diff --git a/cmake/sources-python.cmake b/cmake/sources-python.cmake
index df2ff5a223..5c437f8e42 100644
--- a/cmake/sources-python.cmake
+++ b/cmake/sources-python.cmake
@@ -312,6 +312,7 @@ set(highs_headers_python
     src/lp_data/HighsSolution.h
     src/lp_data/HighsSolutionDebug.h
     src/lp_data/HighsSolve.h
+    src/lp_data/HighsSolverStats.h
     src/lp_data/HighsStatus.h
     src/lp_data/HStruct.h
     src/mip/HighsCliqueTable.h
diff --git a/cmake/sources.cmake b/cmake/sources.cmake
index 96e84b7f5b..4e71dc176c 100644
--- a/cmake/sources.cmake
+++ b/cmake/sources.cmake
@@ -316,6 +316,7 @@ set(highs_headers
     lp_data/HighsSolution.h
     lp_data/HighsSolutionDebug.h
     lp_data/HighsSolve.h
+    lp_data/HighsSolverStats.h
     lp_data/HighsStatus.h
     lp_data/HStruct.h
     mip/HighsCliqueTable.h
diff --git a/src/Highs.h b/src/Highs.h
index 70d800a5de..a382435d8c 100644
--- a/src/Highs.h
+++ b/src/Highs.h
@@ -1216,11 +1216,35 @@ class Highs {
 
   // Start of advanced methods for HiGHS MIP solver
 
-  const HighsSimplexStats& getSimplexStats() const {
-    return ekk_instance_.getSimplexStats();
+  const HighsSimplexStats& getPresolvedLpSimplexStats() const {
+    return presolved_lp_simplex_stats_;
   }
-  void reportSimplexStats(FILE* file) const {
-    ekk_instance_.reportSimplexStats(file);
+
+  void reportPresolvedLpSimplexStats(
+      FILE* file, const HighsInt style = HighsSolverStatsReportPretty) const {
+    presolved_lp_simplex_stats_.report(file, "Presolved LP", style);
+  }
+
+  const HighsSimplexStats& getSimplexStats() const { return simplex_stats_; }
+
+  void reportSimplexStats(
+      FILE* file, const HighsInt style = HighsSolverStatsReportPretty) const {
+    simplex_stats_.report(file, "Original LP", style);
+  }
+
+  void passSimplexStats(const HighsSimplexStats simplex_stats) {
+    this->simplex_stats_ = simplex_stats;
+  }
+
+  const HighsIpxStats& getIpxStats() { return ipx_stats_; }
+
+  void reportIpxStats(FILE* file,
+                      const HighsInt style = HighsSolverStatsReportPretty) {
+    ipx_stats_.report(file, "Original LP", style);
+  }
+
+  void passIpxStats(const HighsIpxStats ipx_stats) {
+    this->ipx_stats_ = ipx_stats;
   }
 
   /**
@@ -1438,6 +1462,9 @@ class Highs {
   HighsSparseMatrix standard_form_matrix_;
 
   HEkk ekk_instance_;
+  HighsSimplexStats simplex_stats_;
+  HighsSimplexStats presolved_lp_simplex_stats_;
+  HighsIpxStats ipx_stats_;
 
   HighsPresolveLog presolve_log_;
 
@@ -1500,7 +1527,8 @@ class Highs {
   //
   // Invalidates all solver data in Highs class members by calling
   // invalidateModelStatus(), invalidateSolution(), invalidateBasis(),
-  // invalidateInfo() and invalidateEkk()
+  // invalidateInfo(), invalidateEkk(), invalidateIis(),
+  // invalidateSimplexStats() and invalidateIpxStats();
   void invalidateUserSolverData();
   //
   // Invalidates the model status, solution_ and info_
@@ -1527,6 +1555,12 @@ class Highs {
   // Invalidates iis_
   void invalidateIis();
 
+  // Invalidates simplex_stats_ and presolved_lp_simplex_stats_;
+  void invalidateSimplexStats();
+
+  // Invalidates ipx_stats_
+  void invalidateIpxStats();
+
   HighsStatus returnFromWriteSolution(FILE* file,
                                       const HighsStatus return_status);
   HighsStatus returnFromRun(const HighsStatus return_status,
diff --git a/src/ipm/IpxWrapper.cpp b/src/ipm/IpxWrapper.cpp
index eae2aff189..44a78ee568 100644
--- a/src/ipm/IpxWrapper.cpp
+++ b/src/ipm/IpxWrapper.cpp
@@ -23,16 +23,17 @@ using std::min;
 
 HighsStatus solveLpIpx(HighsLpSolverObject& solver_object) {
   return solveLpIpx(solver_object.options_, solver_object.timer_,
-                    solver_object.lp_, solver_object.basis_,
-                    solver_object.solution_, solver_object.model_status_,
-                    solver_object.highs_info_, solver_object.callback_);
+                    solver_object.lp_, solver_object.ekk_instance_,
+                    solver_object.basis_, solver_object.solution_,
+                    solver_object.model_status_, solver_object.highs_info_,
+                    solver_object.ipx_stats_, solver_object.callback_);
 }
 
 HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
-                       const HighsLp& lp, HighsBasis& highs_basis,
-                       HighsSolution& highs_solution,
+                       const HighsLp& lp, const HEkk& ekk_instance,
+                       HighsBasis& highs_basis, HighsSolution& highs_solution,
                        HighsModelStatus& model_status, HighsInfo& highs_info,
-                       HighsCallback& callback) {
+                       HighsIpxStats& ipx_stats, HighsCallback& callback) {
   // Use IPX to try to solve the LP
   //
   // Can return HighsModelStatus (HighsStatus) values:
@@ -58,6 +59,7 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   // then a basis and primal+dual solution are obtained.
   //
   //
+  const double entry_run_time = timer.readRunHighsClock();
   // Indicate that there is no valid primal solution, dual solution or basis
   highs_basis.valid = false;
   highs_solution.value_valid = false;
@@ -118,9 +120,13 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   parameters.analyse_basis_data =
       kHighsAnalysisLevelNlaData & options.highs_analysis_level;
   // Determine the run time allowed for IPX
-  parameters.time_limit = options.time_limit - timer.readRunHighsClock();
+  parameters.time_limit = std::max(options.time_limit - entry_run_time, 0.0);
+
   parameters.ipm_maxiter =
       options.ipm_iteration_limit - highs_info.ipm_iteration_count;
+  parameters.cr1_maxiter = options.cr1_iteration_limit;
+  parameters.cr2_maxiter = options.cr2_iteration_limit;
+  parameters.kkt_logging = options.kkt_logging;
   // Determine if crossover is to be run or not
   //
   // When doing analytic centring calculations, crossover must not be
@@ -146,6 +152,8 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   parameters.max_centring_steps = options.max_centring_steps;
   parameters.centring_ratio_tolerance = options.centring_ratio_tolerance;
 
+  parameters.simplex_stats = ekk_instance.getSimplexStats();
+
   // Set the internal IPX parameters
   lps.SetParameters(parameters);
 
@@ -178,11 +186,16 @@ HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
   const bool report_solve_data =
       kHighsAnalysisLevelSolverSummaryData & options.highs_analysis_level;
   // Get solver and solution information.
+  ipx_stats = lps.GetIpxStats();
+  ipx_stats.averages();
+  ipx_stats.valid = true;
   // Struct ipx_info defined in ipx/ipx_info.h
   const ipx::Info ipx_info = lps.GetInfo();
   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 +973,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",
@@ -1074,3 +1091,158 @@ void reportSolveData(const HighsLogOptions& log_options,
               "    Volume increase      = %11.4g\n\n",
               ipx_info.volume_increase);
 }
+
+void HighsIpxStats::workTerms(double* terms) {
+  const double nonbasic_nz =
+      double(this->num_nz + this->num_row - this->average_type2_matrix_nz);
+  terms[HighsIpxWorkTermCr1IterNumRow] = double(this->num_type1_iteration) *
+                                         average_type1_cr_count *
+                                         double(this->num_row);
+  terms[HighsIpxWorkTermCr1IterNumNz] = double(this->num_type1_iteration) *
+                                        average_type1_cr_count *
+                                        double(this->num_nz);
+  terms[HighsIpxWorkTermCr2IterNumRow] = double(this->num_type2_iteration) *
+                                         average_type2_cr_count *
+                                         double(this->num_row);
+  terms[HighsIpxWorkTermCr2IterNumNz] = double(this->num_type2_iteration) *
+                                        average_type2_cr_count *
+                                        (nonbasic_nz + average_type2_invert_nz);
+}
+
+double HighsIpxStats::workEstimate() {
+  double* terms = new double[HighsIpxWorkTermCount];
+  this->workTerms(terms);
+  double work = 0;
+  for (HighsInt iX = 0; iX < HighsIpxWorkTermCount; iX++) {
+    assert(terms[iX] > 0);
+    work += terms[iX] * kIpxWorkCoefficients[iX];
+  }
+  delete[] terms;
+  return work;
+}
+
+void HighsIpxStats::averages() {
+  num_type1_iteration = 0;
+  num_type2_iteration = 0;
+  average_type1_cr_count = 0;
+  average_type2_cr_count = 0;
+  average_type2_matrix_nz = 0;
+  average_type2_invert_nz = 0;
+  for (HighsInt iteration = 0; iteration < this->iteration_count; iteration++) {
+    if (this->cr_type[iteration] == 1) {
+      num_type1_iteration++;
+      average_type1_cr_count += this->cr_count[iteration];
+    } else {
+      num_type2_iteration++;
+      average_type2_cr_count += this->cr_count[iteration];
+      average_type2_matrix_nz += this->factored_basis_num_el[iteration];
+      average_type2_invert_nz += this->invert_num_el[iteration];
+    }
+  }
+  if (num_type1_iteration)
+    average_type1_cr_count /= (1.0 * num_type1_iteration);
+  if (num_type2_iteration) {
+    average_type2_cr_count /= (1.0 * num_type2_iteration);
+    average_type2_matrix_nz /= (1.0 * num_type2_iteration);
+    average_type2_invert_nz /= (1.0 * num_type2_iteration);
+  }
+}
+
+void HighsIpxStats::report(FILE* file, const std::string message,
+                           const HighsInt style) {
+  if (style == HighsSolverStatsReportPretty) {
+    fprintf(file, "\nIpx stats: %s\n", message.c_str());
+    fprintf(file, "   valid                      = %d\n", this->valid);
+    fprintf(file, "   num_col                    = %d\n", this->num_col);
+    fprintf(file, "   num_row                    = %d\n", this->num_row);
+    fprintf(file, "   num_nz                     = %d\n", this->num_nz);
+    fprintf(file, "   iteration_count            = %d\n",
+            this->iteration_count);
+    if (this->iteration_count != HighsInt(cr_type.size()))
+      printf("iteration_count = %d != %d = cr_type.size()\n",
+             int(this->iteration_count), int(this->cr_type.size()));
+    if (this->iteration_count != HighsInt(cr_count.size()))
+      printf("iteration_count = %d != %d = cr_count.size()\n",
+             int(this->iteration_count), int(this->cr_count.size()));
+    if (this->iteration_count != HighsInt(invert_num_el.size()))
+      printf("iteration_count = %d != %d = invert_num_el.size()\n",
+             int(this->iteration_count), int(this->invert_num_el.size()));
+    if (this->iteration_count != HighsInt(factored_basis_num_el.size()))
+      printf("iteration_count = %d != %d = factored_basis_num_el.size()\n",
+             int(this->iteration_count),
+             int(this->factored_basis_num_el.size()));
+    assert(this->iteration_count == HighsInt(this->cr_type.size()));
+    assert(this->iteration_count == HighsInt(this->cr_count.size()));
+    assert(this->iteration_count == HighsInt(this->invert_num_el.size()));
+    assert(this->iteration_count ==
+           HighsInt(this->factored_basis_num_el.size()));
+    if (iteration_count > 0)
+      fprintf(file, "   Iter  type cr_count basisNz invertNz\n");
+    // printf(file, "   dddd     d    ddddd ddddddd  ddddddd\n");
+    for (HighsInt iteration = 0; iteration < iteration_count; iteration++)
+      fprintf(file, "   %4d     %1d    %5d %7d  %7d\n", int(iteration),
+              int(this->cr_type[iteration]), int(this->cr_count[iteration]),
+              int(this->factored_basis_num_el[iteration]),
+              int(this->invert_num_el[iteration]));
+    fprintf(file, "   num_nz                     = %d\n", this->num_nz);
+    fprintf(file, "   Type 1 iteration           = %d\n",
+            this->num_type1_iteration);
+    fprintf(file, "          mean CR count       = %g\n",
+            this->average_type1_cr_count);
+    fprintf(file, "   Type 2 iteration           = %d\n",
+            this->num_type2_iteration);
+    fprintf(file, "          mean CR count       = %g\n",
+            this->average_type2_cr_count);
+    fprintf(file, "          mean CR matrix nz   = %g\n",
+            this->average_type2_matrix_nz);
+    fprintf(file, "          mean CR INVERT nz   = %g\n",
+            this->average_type2_invert_nz);
+    fprintf(file, "   Type 1                time = %g\n", this->type1_time);
+    fprintf(file, "   Starting basis        time = %g\n", this->basis0_time);
+    fprintf(file, "   Type 2                time = %g\n", this->type2_time);
+    fprintf(file, "   IPM                   time = %g\n", this->ipm_time);
+    fprintf(file, "   Crossover             time = %g\n", this->crossover_time);
+  } else if (style == HighsSolverStatsReportCsvHeader) {
+    fprintf(file,
+            "valid,col,row,nz,iteration_count,cr_count,iteration_count, "
+            "cr_count,matrix_nz,invert_nz,type1_time,basis0_time,type2_time,"
+            "ipm_time,crossover_time,");
+  } else if (style == HighsSolverStatsReportCsvData) {
+    this->averages();
+    fprintf(
+        file, "%d,%d,%d,%d,%d,%d,%d,%d,%d,%d,%g,%g,%g,%g,%g,", int(this->valid),
+        int(this->num_col), int(this->num_row), int(this->num_nz),
+        int(this->num_type1_iteration), int(this->average_type1_cr_count),
+        int(this->num_type2_iteration), int(this->average_type2_cr_count),
+        int(this->average_type2_matrix_nz), int(this->average_type2_invert_nz),
+        this->type1_time, this->basis0_time, this->type2_time, this->ipm_time,
+        this->crossover_time);
+  } else {
+    fprintf(file, "Unknown IPX stats report style of %d\n", int(style));
+    assert(123 == 456);
+  }
+}
+
+void HighsIpxStats::initialise() {
+  valid = false;
+  iteration_count = 0;
+  num_col = 0;
+  num_row = 0;
+  num_nz = 0;
+  iteration_count = 0;
+  cr_type.clear();
+  cr_count.clear();
+  factored_basis_num_el.clear();
+  invert_num_el.clear();
+  num_type1_iteration = 0;
+  num_type2_iteration = 0;
+  average_type1_cr_count = 0;
+  average_type2_cr_count = 0;
+  average_type2_matrix_nz = 0;
+  average_type2_invert_nz = 0;
+  type1_time = 0;
+  basis0_time = 0;
+  type2_time = 0;
+  ipm_time = 0;
+  crossover_time = 0;
+}
diff --git a/src/ipm/IpxWrapper.h b/src/ipm/IpxWrapper.h
index 660165d82d..16ccb6df0d 100644
--- a/src/ipm/IpxWrapper.h
+++ b/src/ipm/IpxWrapper.h
@@ -25,10 +25,10 @@
 HighsStatus solveLpIpx(HighsLpSolverObject& solver_object);
 
 HighsStatus solveLpIpx(const HighsOptions& options, HighsTimer& timer,
-                       const HighsLp& lp, HighsBasis& highs_basis,
-                       HighsSolution& highs_solution,
+                       const HighsLp& lp, const HEkk& ekk_instance,
+                       HighsBasis& highs_basis, HighsSolution& highs_solution,
                        HighsModelStatus& model_status, HighsInfo& highs_info,
-                       HighsCallback& callback);
+                       HighsIpxStats& ipx_stats, HighsCallback& callback);
 
 void fillInIpxData(const HighsLp& lp, ipx::Int& num_col, ipx::Int& num_row,
                    std::vector<double>& obj, std::vector<double>& col_lb,
diff --git a/src/ipm/ipx/basiclu_wrapper.cc b/src/ipm/ipx/basiclu_wrapper.cc
index e38686f88c..cc129336bd 100644
--- a/src/ipm/ipx/basiclu_wrapper.cc
+++ b/src/ipm/ipx/basiclu_wrapper.cc
@@ -253,6 +253,14 @@ bool BasicLu::_NeedFreshFactorization() {
     return nforrest == dim || update_cost > 1.0;
 }
 
+Int BasicLu::_matrix_nz() const {
+  return xstore_[BASICLU_MATRIX_NZ];
+}
+
+Int BasicLu::_invert_nz() const {
+  return xstore_[BASICLU_LNZ] + xstore_[BASICLU_UNZ] + xstore_[BASICLU_DIM];
+}
+
 double BasicLu::_fill_factor() const {
     return fill_factor_;
 }
diff --git a/src/ipm/ipx/basiclu_wrapper.h b/src/ipm/ipx/basiclu_wrapper.h
index 6ad9c5775d..b2e04a903b 100644
--- a/src/ipm/ipx/basiclu_wrapper.h
+++ b/src/ipm/ipx/basiclu_wrapper.h
@@ -24,6 +24,8 @@ class BasicLu : public LuUpdate {
     void _BtranForUpdate(Int j, IndexedVector& lhs) override;
     Int _Update(double pivot) override;
     bool _NeedFreshFactorization() override;
+    Int _matrix_nz() const override;
+    Int _invert_nz() const override;
     double _fill_factor() const override;
     double _pivottol() const override;
     void _pivottol(double new_pivottol) override;
@@ -39,6 +41,8 @@ class BasicLu : public LuUpdate {
     std::vector<double> xstore_;
     std::vector<Int> Li_, Ui_, Wi_;
     std::vector<double> Lx_, Ux_, Wx_;
+    Int matrix_nz_;
+    Int invert_nz_;
     double fill_factor_;
 };
 
diff --git a/src/ipm/ipx/basis.cc b/src/ipm/ipx/basis.cc
index 607de407d1..0275d9945b 100644
--- a/src/ipm/ipx/basis.cc
+++ b/src/ipm/ipx/basis.cc
@@ -130,6 +130,9 @@ Int Basis::Factorize() {
         Int flag = lu_->Factorize(begin.data(), end.data(), AI.rowidx(),
                                   AI.values(), false);
         num_factorizations_++;
+	matrix_nz_ = lu_->matrix_nz();
+	invert_nz_ = lu_->invert_nz();
+	double fill_factor = lu_->fill_factor();
         fill_factors_.push_back(lu_->fill_factor());
         if (flag & 2) {
             AdaptToSingularFactorization();
@@ -449,6 +452,21 @@ double Basis::time_update() const {
     return time_update_;
 }
 
+Int Basis::matrix_nz() const {
+  return matrix_nz_;
+}
+
+Int Basis::invert_nz() const {
+  return invert_nz_;
+}
+
+double Basis::current_fill() const {
+    if (fill_factors_.empty())
+        return 0.0;
+    Int num_factors = fill_factors_.size();
+    return fill_factors_[num_factors-1];
+}
+
 double Basis::mean_fill() const {
     if (fill_factors_.empty())
         return 0.0;
@@ -630,6 +648,8 @@ void Basis::CrashFactorize(Int* num_dropped) {
     Int flag = lu_->Factorize(begin.data(), end.data(), AI.rowidx(),
                               AI.values(), true);
     num_factorizations_++;
+    matrix_nz_ = lu_->matrix_nz();
+    invert_nz_ = lu_->invert_nz();
     fill_factors_.push_back(lu_->fill_factor());
     Int ndropped = 0;
     if (flag & 2)
diff --git a/src/ipm/ipx/basis.h b/src/ipm/ipx/basis.h
index 4ae358c171..82cdd27749 100644
--- a/src/ipm/ipx/basis.h
+++ b/src/ipm/ipx/basis.h
@@ -219,6 +219,9 @@ class Basis {
     double time_ftran() const;        // time FTRAN, including partial
     double time_btran() const;        // time BTRAN, including partial
     double time_update() const;       // time LU update
+    Int matrix_nz() const;            // Current nonzeros of matrix to be factored
+    Int invert_nz() const;            // Nonzeros of factored matrix
+    double current_fill() const;      // Current LU fill factors
     double mean_fill() const;         // geom. mean of LU fill factors
     double max_fill() const;          // max LU fill factor
     
@@ -307,6 +310,8 @@ class Basis {
     double time_btran_{0.0};       // time for BTRAN ops, including partial
     double time_update_{0.0};      // time for LU updates
     double time_factorize_{0.0};   // time for LU factorizations
+    Int matrix_nz_{0};             // nonzeros of matrix to be factored
+    Int invert_nz_{0};             // nonzeros of INVERT of factored matrix
     std::vector<double> fill_factors_; // fill factors from LU factorizations
     double sum_ftran_density_{0.0};
     double sum_btran_density_{0.0};
diff --git a/src/ipm/ipx/conjugate_residuals.cc b/src/ipm/ipx/conjugate_residuals.cc
index 9fe770e3b9..26b89411b1 100644
--- a/src/ipm/ipx/conjugate_residuals.cc
+++ b/src/ipm/ipx/conjugate_residuals.cc
@@ -12,6 +12,9 @@ ConjugateResiduals::ConjugateResiduals(const Control& control) :
 void ConjugateResiduals::Solve(LinearOperator& C, const Vector& rhs,
                               double tol, const double* resscale, Int maxiter,
                               Vector& lhs) {
+
+    const bool cr_logging = control_.kkt_logging();
+
     const Int m = rhs.size();
     Vector residual(m);  // rhs - C*lhs
     Vector step(m);      // update to lhs
@@ -23,9 +26,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 +43,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 +55,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 +67,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,15 +86,18 @@ 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();
 }
 
 void ConjugateResiduals::Solve(LinearOperator& C, LinearOperator& P,
                               const Vector& rhs, double tol,
                               const double* resscale, Int maxiter, Vector& lhs){
+    const bool cr_logging = control_.kkt_logging();
+
     const Int m = rhs.size();
     Vector residual(m);   // rhs - C*lhs
     Vector sresidual(m);  // preconditioned residual
@@ -109,9 +119,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 +137,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 +149,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 +222,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/ipm/ipx/control.cc b/src/ipm/ipx/control.cc
index 9a149fab15..4b6a3c58e5 100644
--- a/src/ipm/ipx/control.cc
+++ b/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_);
diff --git a/src/ipm/ipx/control.h b/src/ipm/ipx/control.h
index d322b0d383..d70933eae9 100644
--- a/src/ipm/ipx/control.h
+++ b/src/ipm/ipx/control.h
@@ -68,7 +68,10 @@ class Control {
     double ipm_optimality_tol() const { return parameters_.ipm_optimality_tol; }
     double ipm_drop_primal() const { return parameters_.ipm_drop_primal; }
     double ipm_drop_dual() const { return parameters_.ipm_drop_dual; }
+    bool kkt_logging() const { return parameters_.kkt_logging; }
     double kkt_tol() const { return parameters_.kkt_tol; }
+    ipxint cr1_maxiter() const { return parameters_.cr1_maxiter; }
+    ipxint cr2_maxiter() const { return parameters_.cr2_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; }
@@ -91,8 +94,10 @@ class Control {
     double centringRatioReduction() const {return parameters_.centring_ratio_reduction; }
     double centringAlphaScaling() const{return parameters_.centring_alpha_scaling; }
     ipxint badProductsTolerance() const{return parameters_.bad_products_tolerance; }
+    HighsSimplexStats simplexStats() const{return parameters_.simplex_stats; }
+    HighsIpxStats ipxStats() const{return parameters_.ipx_stats; }
 
-    const Parameters& parameters() const;
+   const Parameters& parameters() const;
     void parameters(const Parameters& new_parameters);
     void callback(HighsCallback* callback);
 
diff --git a/src/ipm/ipx/forrest_tomlin.cc b/src/ipm/ipx/forrest_tomlin.cc
index 96004f8664..4d7f38c7a1 100644
--- a/src/ipm/ipx/forrest_tomlin.cc
+++ b/src/ipm/ipx/forrest_tomlin.cc
@@ -235,6 +235,10 @@ bool ForrestTomlin::_NeedFreshFactorization() {
     return false;
 }
 
+Int ForrestTomlin::_matrix_nz() const { return 0;}
+
+Int ForrestTomlin::_invert_nz() const { return 0;}
+
 double ForrestTomlin::_fill_factor() const {
     return fill_factor_;
 }
diff --git a/src/ipm/ipx/forrest_tomlin.h b/src/ipm/ipx/forrest_tomlin.h
index dc60bfd6ab..73a3134daf 100644
--- a/src/ipm/ipx/forrest_tomlin.h
+++ b/src/ipm/ipx/forrest_tomlin.h
@@ -44,6 +44,8 @@ class ForrestTomlin : public LuUpdate {
     void _BtranForUpdate(Int j, IndexedVector& lhs) override;
     Int _Update(double pivot) override;
     bool _NeedFreshFactorization() override;
+    Int _matrix_nz() const override;
+    Int _invert_nz() const override;
     double _fill_factor() const override;
     double _pivottol() const override;
     void _pivottol(double new_pivottol) override;
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/ipm.cc b/src/ipm/ipx/ipm.cc
index 4a21364951..fca4730476 100644
--- a/src/ipm/ipx/ipm.cc
+++ b/src/ipm/ipx/ipm.cc
@@ -20,10 +20,11 @@ struct IPM::Step {
 
 IPM::IPM(const Control& control) : control_(control) {}
 
-void IPM::StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info) {
+  void IPM::StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info, HighsIpxStats* ipx_stats) {
     kkt_ = kkt;
     iterate_ = iterate;
     info_ = info;
+    ipx_stats_ = ipx_stats;
     PrintHeader();
     ComputeStartingPoint();
     if (info->errflag == 0)
@@ -42,7 +43,7 @@ void IPM::StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info) {
     }
 }
 
-void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
+void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info, const bool diag) {
     const Model& model = iterate->model();
     const Int m = model.rows();
     const Int n = model.cols();
@@ -86,6 +87,7 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
         }
         if ((info->errflag = control_.InterruptCheck(info->iter)) != 0)
             break;
+	Int kktiter = diag ? -info_->kktiter1 : -info_->kktiter2;
         kkt->Factorize(iterate, info);
         if (info->errflag)
             break;
@@ -96,7 +98,23 @@ void IPM::Driver(KKTSolver* kkt, Iterate* iterate, Info* info) {
         if (info->errflag)
             break;
         MakeStep(step);
+	//	HighsSimplexStats simplex_stats = control_.simplexStats();
+	//	double simplex_work_measure = simplex_stats.workEstimate();
+	//
         info->iter++;
+	
+	// Update IPX stats
+	kktiter += diag ? info_->kktiter1 : info_->kktiter2;
+	Int cr_type = diag ? 1 : 2;
+	ipx_stats_->iteration_count++;
+	ipx_stats_->cr_type.push_back(cr_type);
+	ipx_stats_->cr_count.push_back(kktiter);
+	Int matrix_nz = diag ? 0 : kkt_->basis()->matrix_nz();
+	Int invert_nz = diag ? 0 : kkt_->basis()->invert_nz();
+	ipx_stats_->factored_basis_num_el.push_back(matrix_nz);
+	ipx_stats_->invert_num_el.push_back(invert_nz);
+	//	if (cr_type == 2) ipx_stats_->report(stdout, "IPM::Driver()");
+
         PrintOutput();
     }
 
@@ -202,7 +220,9 @@ void IPM::ComputeStartingPoint() {
     Vector x(n+m), xl(n+m), xu(n+m), y(m), zl(n+m), zu(n+m);
     Vector rb(m);               // workspace
 
-    // Factorize the KKT matrix with the identity matrix in the (1,1) block.
+    // "Factorize" the KKT matrix with the identity matrix in the (1,1) block.
+    //
+    // Just inverts the diagonal of the normal matrix
     kkt_->Factorize(nullptr, info_);
     if (info_->errflag)
         return;
@@ -222,6 +242,7 @@ void IPM::ComputeStartingPoint() {
     }
     double tol = 0.1 * Infnorm(rb);
     zl = 0.0;
+    Int kktiter1 = -info_->kktiter1;
     kkt_->Solve(zl, rb, tol, xl, y, info_);
     if (info_->errflag)
         return;
@@ -324,6 +345,15 @@ void IPM::ComputeStartingPoint() {
     }
     iterate_->Initialize(x, xl, xu, y, zl, zu);
     best_complementarity_ = iterate_->complementarity();
+
+    // Update IPX stats
+    kktiter1 += info_->kktiter1;
+    ipx_stats_->iteration_count++;
+    ipx_stats_->cr_type.push_back(1);
+    ipx_stats_->cr_count.push_back(kktiter1);
+    ipx_stats_->factored_basis_num_el.push_back(0);
+    ipx_stats_->invert_num_el.push_back(0);
+    //    ipx_stats_->report(stdout, "IPM::ComputeStartingPoint()");
 }
 
 // Computes maximum alpha such that x + alpha*dx >= 0.
@@ -856,7 +886,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);
diff --git a/src/ipm/ipx/ipm.h b/src/ipm/ipx/ipm.h
index 7d5ae83e60..ab5f2786b9 100644
--- a/src/ipm/ipx/ipm.h
+++ b/src/ipm/ipx/ipm.h
@@ -25,7 +25,7 @@ class IPM {
     // IPX_STATUS_time_limit if the KKT solver was interrupted by time limit,
     // IPX_STATUS_failed     if the KKT solver failed with info->errflag.
     // If the method did not terminate successfully, @iterate is unchanged.
-    void StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info);
+    void StartingPoint(KKTSolver* kkt, Iterate* iterate, Info* info, HighsIpxStats* ipx_stats);
 
     // Updates @iterate by interior point iterations. On return ipm_status is
     // IPX_STATUS_optimal       if iterate->term_crit_reached() is true,
@@ -33,7 +33,7 @@ class IPM {
     // IPX_STATUS_no_progress   if no progress over a number of iterations,
     // IPX_STATUS_time_limit    if interrupted by time limit,
     // IPX_STATUS_failed        if the KKT solver failed with info->errflag.
-    void Driver(KKTSolver* kkt, Iterate* iterate, Info* info);
+    void Driver(KKTSolver* kkt, Iterate* iterate, Info* info, const bool diag);
 
     Int maxiter() const { return maxiter_; }
     void maxiter(Int i) { maxiter_ = i; }
@@ -74,6 +74,7 @@ class IPM {
     KKTSolver* kkt_{nullptr};
     Iterate* iterate_{nullptr};
     Info* info_{nullptr};
+    HighsIpxStats* ipx_stats_{nullptr};
 
     double step_primal_{0.0}, step_dual_{0.0};
     // Counts the # bad iterations since the last good iteration. An iteration
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/ipx_internal.h b/src/ipm/ipx/ipx_internal.h
index f51654b2ac..4fe05ea86f 100644
--- a/src/ipm/ipx/ipx_internal.h
+++ b/src/ipm/ipx/ipx_internal.h
@@ -31,7 +31,10 @@ struct Parameters : public ipx_parameters {
     ipm_optimality_tol = 1e-8;
     ipm_drop_primal = 1e-9;
     ipm_drop_dual = 1e-9;
+    kkt_logging = false;
     kkt_tol = 0.3;
+    cr1_maxiter = -1;
+    cr2_maxiter = -1;
     crash_basis = 1;
     dependency_tol = 1e-6;
     volume_tol = 2.0;
diff --git a/src/ipm/ipx/ipx_parameters.h b/src/ipm/ipx/ipx_parameters.h
index 2f9db89186..5ea11bc3aa 100644
--- a/src/ipm/ipx/ipx_parameters.h
+++ b/src/ipm/ipx/ipx_parameters.h
@@ -2,6 +2,7 @@
 #define IPX_PARAMETERS_H_
 
 #include "io/HighsIO.h"
+#include "lp_data/HighsSolverStats.h"
 #include "ipm/ipx/ipx_config.h"
 #include <iostream>
 
@@ -28,6 +29,9 @@ struct ipx_parameters {
     double ipm_drop_dual;
 
     /* Linear solver */
+    bool kkt_logging;
+    ipxint cr1_maxiter;
+    ipxint cr2_maxiter;
     double kkt_tol;
 
     /* Basis construction in IPM */
@@ -65,7 +69,10 @@ struct ipx_parameters {
     /* HiGHS logging parameters */
     bool highs_logging;
     const HighsLogOptions* log_options;
-  
+
+    /* Simplex and IPX solution stats */
+    HighsSimplexStats simplex_stats;
+    HighsIpxStats ipx_stats;
 };
 
 #ifdef __cplusplus
diff --git a/src/ipm/ipx/kkt_solver.cc b/src/ipm/ipx/kkt_solver.cc
index dde0967e52..a828fcb0a3 100644
--- a/src/ipm/ipx/kkt_solver.cc
+++ b/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(); }
 
diff --git a/src/ipm/ipx/kkt_solver.h b/src/ipm/ipx/kkt_solver.h
index 3889f1b97c..e5ab15bc22 100644
--- a/src/ipm/ipx/kkt_solver.h
+++ b/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; }
 };
diff --git a/src/ipm/ipx/kkt_solver_basis.cc b/src/ipm/ipx/kkt_solver_basis.cc
index dddb35ec8b..46d28e2766 100644
--- a/src/ipm/ipx/kkt_solver_basis.cc
+++ b/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++)
diff --git a/src/ipm/ipx/kkt_solver_basis.h b/src/ipm/ipx/kkt_solver_basis.h
index db7e1af131..e2cf452de6 100644
--- a/src/ipm/ipx/kkt_solver_basis.h
+++ b/src/ipm/ipx/kkt_solver_basis.h
@@ -33,7 +33,9 @@ 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 +59,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};
 };
 
diff --git a/src/ipm/ipx/kkt_solver_diag.cc b/src/ipm/ipx/kkt_solver_diag.cc
index 382c9d06cc..9e66ad6848 100644
--- a/src/ipm/ipx/kkt_solver_diag.cc
+++ b/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++)
diff --git a/src/ipm/ipx/kkt_solver_diag.h b/src/ipm/ipx/kkt_solver_diag.h
index 31a2b708ba..0b9ae559ae 100644
--- a/src/ipm/ipx/kkt_solver_diag.h
+++ b/src/ipm/ipx/kkt_solver_diag.h
@@ -24,12 +24,13 @@ class KKTSolverDiag : public KKTSolver {
 
     Int maxiter() const { return maxiter_; }
     void maxiter(Int new_maxiter) { maxiter_ = new_maxiter; }
-
+  
 private:
     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
diff --git a/src/ipm/ipx/lp_solver.cc b/src/ipm/ipx/lp_solver.cc
index 145930c75b..b3ec9b6f3d 100644
--- a/src/ipm/ipx/lp_solver.cc
+++ b/src/ipm/ipx/lp_solver.cc
@@ -55,7 +55,9 @@ Int LpSolver::Solve() {
     control_.OpenLogfile();
     control_.hLog("IPX version 1.0\n");
     try {
+        ipx_stats_.ipm_time = -control_.Elapsed();
         InteriorPointSolve();
+        ipx_stats_.ipm_time += control_.Elapsed();
 	const bool run_crossover_on = control_.run_crossover() == 1;
 	const bool run_crossover_choose = control_.run_crossover() == -1;
 	const bool run_crossover_not_off = run_crossover_choose || run_crossover_on;
@@ -65,16 +67,18 @@ 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);
+	  }
+	  ipx_stats_.crossover_time = -control_.Elapsed();
+	  BuildCrossoverStartingPoint();
+	  RunCrossover();
+	  ipx_stats_.crossover_time += control_.Elapsed();
         }
         if (basis_) {
             info_.ftran_sparse = basis_->frac_ftran_sparse();
@@ -128,6 +132,10 @@ Info LpSolver::GetInfo() const {
     return info_;
 }
 
+HighsIpxStats LpSolver::GetIpxStats() const {
+    return ipx_stats_;
+}
+
 Int LpSolver::GetInteriorSolution(double* x, double* xl, double* xu,
                                   double* slack, double* y, double* zl,
                                   double* zu) const {
@@ -392,24 +400,37 @@ void LpSolver::RunIPM() {
     info_.centring_tried = false;
     info_.centring_success = false;
 
+    ipx_stats_.initialise();
+    ipx_stats_.num_col = model_.cols();
+    ipx_stats_.num_row = model_.rows();
+    ipx_stats_.num_nz = model_.AI().entries();
+    
     if (x_start_.size() != 0) {
       control_.hLog(" Using starting point provided by user. Skipping initial iterations.\n");
         iterate_->Initialize(x_start_, xl_start_, xu_start_,
                              y_start_, zl_start_, zu_start_);
     }
     else {
+        ipx_stats_.type1_time = -control_.Elapsed();
         ComputeStartingPoint(ipm);
-        if (info_.status_ipm != IPX_STATUS_not_run)
-            return;
+        if (info_.status_ipm != IPX_STATUS_not_run) {
+	  ipx_stats_.type1_time += control_.Elapsed();
+	  return;
+	}
         RunInitialIPM(ipm);
+        ipx_stats_.type1_time += control_.Elapsed();
         if (info_.status_ipm != IPX_STATUS_not_run)
             return;
     }
+    ipx_stats_.basis0_time = -control_.Elapsed();
     BuildStartingBasis();
+    ipx_stats_.basis0_time += control_.Elapsed();
     if (info_.status_ipm != IPX_STATUS_not_run || 
         info_.centring_tried)
         return;
+    ipx_stats_.type2_time = -control_.Elapsed();
     RunMainIPM(ipm);
+    ipx_stats_.type2_time += control_.Elapsed();
 }
 
 void LpSolver::MakeIPMStartingPointValid() {
@@ -474,7 +495,7 @@ void LpSolver::ComputeStartingPoint(IPM& ipm) {
     // If the starting point procedure fails, then iterate_ remains as
     // initialized by the constructor, which is a valid state for
     // postprocessing/postsolving.
-    ipm.StartingPoint(&kkt, iterate_.get(), &info_);
+    ipm.StartingPoint(&kkt, iterate_.get(), &info_, &ipx_stats_);
     info_.time_ipm1 += timer.Elapsed();
 }
 
@@ -484,15 +505,28 @@ void LpSolver::RunInitialIPM(IPM& ipm) {
 
     Int switchiter = control_.switchiter();
     if (switchiter < 0) {
-        // Switch iteration not specified by user. Run as long as KKT solver
-        // converges within min(500,10+m/20) iterations.
+        // Switch iteration not specified by user. Run as long as KKT
+        // solver converges within
+        //
+        // min(control_.cr1_maxiter(), 500,10+m/20)
+        //
+        // iterations, ignoring control_.cr1_maxiter() if negative
+      //
+      // 2049 revert all this
         Int m = model_.rows();
-        kkt.maxiter(std::min(500l, (long) (10+m/20) ));
+	ipxint df_cr1_maxiter = std::min(500l, (long) (10+m/20));
+	ipxint cr1_maxiter = df_cr1_maxiter;
+	ipxint control_cr1_maxiter = control_.cr1_maxiter();
+	if (control_cr1_maxiter > 0) cr1_maxiter = std::min(control_cr1_maxiter, cr1_maxiter);
+	//	printf("LpSolver::RunInitialIPM Using kkt.maxiter = %d, from df_cr1_maxiter = %d and control_cr1_maxiter = %d\n",
+	//	       int(cr1_maxiter), int(df_cr1_maxiter), int(control_cr1_maxiter));
+        kkt.maxiter(cr1_maxiter);
         ipm.maxiter(control_.ipm_maxiter());
     } else {
         ipm.maxiter(std::min(switchiter, control_.ipm_maxiter()));
     }
-    ipm.Driver(&kkt, iterate_.get(), &info_);
+    const bool diag = true;
+    ipm.Driver(&kkt, iterate_.get(), &info_, diag);
     switch (info_.status_ipm) {
     case IPX_STATUS_optimal:
         // If the IPM reached its termination criterion in the initial
@@ -558,7 +592,9 @@ void LpSolver::RunMainIPM(IPM& ipm) {
     KKTSolverBasis kkt(control_, *basis_);
     Timer timer;
     ipm.maxiter(control_.ipm_maxiter());
-    ipm.Driver(&kkt, iterate_.get(), &info_);
+    kkt.maxiter(control_.cr2_maxiter());
+    const bool diag = false;
+    ipm.Driver(&kkt, iterate_.get(), &info_, diag);
     info_.time_ipm2 = timer.Elapsed();
 }
 
diff --git a/src/ipm/ipx/lp_solver.h b/src/ipm/ipx/lp_solver.h
index 18eb70606f..1c97ffb27e 100644
--- a/src/ipm/ipx/lp_solver.h
+++ b/src/ipm/ipx/lp_solver.h
@@ -75,6 +75,9 @@ class LpSolver {
     // documentation for the meaning of Info values.
     Info GetInfo() const;
 
+    // Returns the IPX stats 
+    HighsIpxStats GetIpxStats() const;
+
     // Returns the final IPM iterate from the last call to Solve() into user
     // arrays. An iterate is available if GetInfo().status_ipm !=
     // IPX_STATUS_not_run. If no iterate is available, the method does nothing.
@@ -163,6 +166,8 @@ class LpSolver {
     // Returns -1 if no basis was available and 0 otherwise.
     Int SymbolicInvert(Int* rowcounts, Int* colcounts);
 
+    HighsIpxStats ipx_stats_;
+
 private:
     void ClearSolution();
     void InteriorPointSolve();
diff --git a/src/ipm/ipx/lu_update.cc b/src/ipm/ipx/lu_update.cc
index bd2d264319..a1734dc129 100644
--- a/src/ipm/ipx/lu_update.cc
+++ b/src/ipm/ipx/lu_update.cc
@@ -43,6 +43,10 @@ bool LuUpdate::NeedFreshFactorization() {
     return _NeedFreshFactorization();
 }
 
+Int LuUpdate::matrix_nz() const { return _matrix_nz(); }
+
+Int LuUpdate::invert_nz() const { return _invert_nz(); }
+
 double LuUpdate::fill_factor() const {
     return _fill_factor();
 }
diff --git a/src/ipm/ipx/lu_update.h b/src/ipm/ipx/lu_update.h
index 299ec4d8c0..87cd41b6a9 100644
--- a/src/ipm/ipx/lu_update.h
+++ b/src/ipm/ipx/lu_update.h
@@ -94,6 +94,12 @@ class LuUpdate {
     // or speed. No stability test is done.
     bool NeedFreshFactorization();
 
+    // Returns nnz(B) from the last factorization.
+    Int matrix_nz() const;
+
+    // Returns (nnz(L)+nnz(U)) from the last factorization.
+    Int invert_nz() const;
+
     // Returns (nnz(L)+nnz(U))/nnz(B) from the last factorization.
     double fill_factor() const;
 
@@ -117,6 +123,8 @@ class LuUpdate {
     virtual void _BtranForUpdate(Int p, IndexedVector& lhs) = 0;
     virtual Int _Update(double pivot) = 0;
     virtual bool _NeedFreshFactorization() = 0;
+    virtual Int _matrix_nz() const = 0;
+    virtual Int _invert_nz() const = 0;
     virtual double _fill_factor() const = 0;
     virtual double _pivottol() const = 0;
     virtual void _pivottol(double new_pivottol) = 0;
diff --git a/src/lp_data/HConst.h b/src/lp_data/HConst.h
index 5e6a9d67fc..fcea91b853 100644
--- a/src/lp_data/HConst.h
+++ b/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 = -1,
+  kMipAnalyticCentreCalulationChoose = kMipAnalyticCentreCalulationMin,
+  kMipAnalyticCentreCalulationOff,
+  kMipAnalyticCentreCalulationOn,
+  kMipAnalyticCentreCalulationMax = kMipAnalyticCentreCalulationOn
+};
+
 enum class HighsLogType { kInfo = 1, kDetailed, kVerbose, kWarning, kError };
 
 enum SimplexScaleStrategy {
diff --git a/src/lp_data/HStruct.h b/src/lp_data/HStruct.h
index 35a0b0ee42..3963d9df1c 100644
--- a/src/lp_data/HStruct.h
+++ b/src/lp_data/HStruct.h
@@ -154,18 +154,4 @@ struct HighsLinearObjective {
   void clear();
 };
 
-struct HighsSimplexStats {
-  bool valid;
-  HighsInt iteration_count;
-  HighsInt num_invert;
-  HighsInt last_invert_num_el;
-  HighsInt last_factored_basis_num_el;
-  double col_aq_density;
-  double row_ep_density;
-  double row_ap_density;
-  double row_DSE_density;
-  void report(FILE* file, const std::string message = "") const;
-  void initialise(const HighsInt iteration_count_ = 0);
-};
-
 #endif /* LP_DATA_HSTRUCT_H_ */
diff --git a/src/lp_data/Highs.cpp b/src/lp_data/Highs.cpp
index 6352a68a09..c2c18619af 100644
--- a/src/lp_data/Highs.cpp
+++ b/src/lp_data/Highs.cpp
@@ -1275,6 +1275,8 @@ HighsStatus Highs::solve() {
     solveLp(incumbent_lp,
             "Solving LP without presolve, or with basis, or unconstrained",
             this_solve_original_lp_time);
+    simplex_stats_ = this->ekk_instance_.getSimplexStats();
+    simplex_stats_.simplex_time = this_solve_original_lp_time;
     return_status = interpretCallStatus(options_.log_options, call_status,
                                         return_status, "callSolveLp");
     if (return_status == HighsStatus::kError)
@@ -1318,6 +1320,9 @@ HighsStatus Highs::solve() {
         ekk_instance_.lp_name_ = "Original LP";
         solveLp(incumbent_lp, "Not presolved: solving the LP",
                 this_solve_original_lp_time);
+        simplex_stats_ = this->ekk_instance_.getSimplexStats();
+        simplex_stats_.simplex_time = this_solve_original_lp_time;
+        presolved_lp_simplex_stats_ = simplex_stats_;
         return_status = interpretCallStatus(options_.log_options, call_status,
                                             return_status, "callSolveLp");
         if (return_status == HighsStatus::kError)
@@ -1330,6 +1335,9 @@ HighsStatus Highs::solve() {
         reportPresolveReductions(log_options, incumbent_lp, false);
         solveLp(incumbent_lp, "Problem not reduced by presolve: solving the LP",
                 this_solve_original_lp_time);
+        simplex_stats_ = this->ekk_instance_.getSimplexStats();
+        simplex_stats_.simplex_time = this_solve_original_lp_time;
+        presolved_lp_simplex_stats_ = simplex_stats_;
         return_status = interpretCallStatus(options_.log_options, call_status,
                                             return_status, "callSolveLp");
         if (return_status == HighsStatus::kError)
@@ -1378,6 +1386,8 @@ HighsStatus Highs::solve() {
         options_.objective_bound = kHighsInf;
         solveLp(reduced_lp, "Solving the presolved LP",
                 this_solve_presolved_lp_time);
+        presolved_lp_simplex_stats_ = this->ekk_instance_.getSimplexStats();
+        presolved_lp_simplex_stats_.simplex_time = this_solve_presolved_lp_time;
         if (ekk_instance_.status_.initialised_for_solve) {
           // Record the pivot threshold resulting from solving the presolved LP
           // with simplex
@@ -1442,6 +1452,8 @@ HighsStatus Highs::solve() {
                 "Solving the original LP with primal simplex "
                 "to determine infeasible or unbounded",
                 this_solve_original_lp_time);
+        simplex_stats_ = this->ekk_instance_.getSimplexStats();
+        simplex_stats_.simplex_time = this_solve_original_lp_time;
         // Recover the options
         options_ = save_options;
         if (return_status == HighsStatus::kError)
@@ -1579,6 +1591,8 @@ HighsStatus Highs::solve() {
             solveLp(incumbent_lp,
                     "Solving the original LP from the solution after postsolve",
                     this_solve_original_lp_time);
+            simplex_stats_ = this->ekk_instance_.getSimplexStats();
+            simplex_stats_.simplex_time = this_solve_original_lp_time;
             // Determine the iteration count
             postsolve_iteration_count += info_.simplex_iteration_count;
             return_status =
@@ -2326,8 +2340,8 @@ HighsStatus Highs::setBasis(const HighsBasis& basis,
       HighsBasis modifiable_basis = basis;
       modifiable_basis.was_alien = true;
       HighsLpSolverObject solver_object(model_.lp_, modifiable_basis, solution_,
-                                        info_, ekk_instance_, callback_,
-                                        options_, timer_);
+                                        info_, ekk_instance_, ipx_stats_,
+                                        callback_, options_, timer_);
       HighsStatus return_status = formSimplexLpBasisAndFactor(solver_object);
       if (return_status != HighsStatus::kOk) return HighsStatus::kError;
       // Update the HiGHS basis
@@ -3583,6 +3597,8 @@ void Highs::invalidateUserSolverData() {
   invalidateInfo();
   invalidateEkk();
   invalidateIis();
+  invalidateSimplexStats();
+  invalidateIpxStats();
 }
 
 void Highs::invalidateModelStatusSolutionAndInfo() {
@@ -3622,6 +3638,13 @@ void Highs::invalidateEkk() { ekk_instance_.invalidate(); }
 
 void Highs::invalidateIis() { iis_.invalidate(); }
 
+void Highs::invalidateSimplexStats() {
+  simplex_stats_.initialise();
+  presolved_lp_simplex_stats_.initialise();
+}
+
+void Highs::invalidateIpxStats() { ipx_stats_.initialise(); }
+
 HighsStatus Highs::completeSolutionFromDiscreteAssignment() {
   // Determine whether the current solution of a MIP is feasible and,
   // if not, try to assign values to continuous variables and discrete
@@ -3764,7 +3787,7 @@ HighsStatus Highs::callSolveLp(HighsLp& lp, const string message) {
   HighsStatus return_status = HighsStatus::kOk;
 
   HighsLpSolverObject solver_object(lp, basis_, solution_, info_, ekk_instance_,
-                                    callback_, options_, timer_);
+                                    ipx_stats_, callback_, options_, timer_);
 
   // Check that the model is column-wise
   assert(model_.lp_.a_matrix_.isColwise());
@@ -4191,6 +4214,7 @@ HighsStatus Highs::callRunPostsolve(const HighsSolution& solution,
             "Solving the original LP from the solution after postsolve");
         // Determine the timing record
         timer_.stop(timer_.solve_clock);
+        simplex_stats_ = this->ekk_instance_.getSimplexStats();
         return_status = interpretCallStatus(options_.log_options, call_status,
                                             return_status, "callSolveLp");
         // Recover the options
diff --git a/src/lp_data/HighsInfo.cpp b/src/lp_data/HighsInfo.cpp
index 2792f64550..bbb07eb170 100644
--- a/src/lp_data/HighsInfo.cpp
+++ b/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) {
diff --git a/src/lp_data/HighsInfo.h b/src/lp_data/HighsInfo.h
index 3292b8281e..70b56c20ef 100644
--- a/src/lp_data/HighsInfo.h
+++ b/src/lp_data/HighsInfo.h
@@ -116,6 +116,8 @@ struct HighsInfoStruct {
   double max_complementarity_violation;
   double sum_complementarity_violations;
   double primal_dual_integral;
+  HighsInt max_cr_iteration_count1;
+  HighsInt max_cr_iteration_count2;
 };
 
 class HighsInfo : public HighsInfoStruct {
@@ -279,6 +281,18 @@ class HighsInfo : public HighsInfoStruct {
         new InfoRecordDouble("primal_dual_integral", "Primal-dual integral",
                              advanced, &primal_dual_integral, 0);
     records.push_back(record_double);
+
+    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);
   }
 
  public:
diff --git a/src/lp_data/HighsInterface.cpp b/src/lp_data/HighsInterface.cpp
index 1c7ad16b9f..126dbe02f9 100644
--- a/src/lp_data/HighsInterface.cpp
+++ b/src/lp_data/HighsInterface.cpp
@@ -1435,8 +1435,8 @@ HighsStatus Highs::getBasicVariablesInterface(HighsInt* basic_variables) {
     // The LP has no invert to use, so have to set one up, but only
     // for the current basis, so return_value is the rank deficiency.
     HighsLpSolverObject solver_object(lp, basis_, solution_, info_,
-                                      ekk_instance_, callback_, options_,
-                                      timer_);
+                                      ekk_instance_, ipx_stats_, callback_,
+                                      options_, timer_);
     const bool only_from_known_basis = true;
     return_status = interpretCallStatus(
         options_.log_options,
@@ -1945,7 +1945,8 @@ HighsStatus Highs::getPrimalRayInterface(bool& has_primal_ray,
 
 HighsStatus Highs::getRangingInterface() {
   HighsLpSolverObject solver_object(model_.lp_, basis_, solution_, info_,
-                                    ekk_instance_, callback_, options_, timer_);
+                                    ekk_instance_, ipx_stats_, callback_,
+                                    options_, timer_);
   solver_object.model_status_ = model_status_;
   return getRangingData(this->ranging_, solver_object);
 }
diff --git a/src/lp_data/HighsLpSolverObject.h b/src/lp_data/HighsLpSolverObject.h
index 895922487a..d47ca65e0f 100644
--- a/src/lp_data/HighsLpSolverObject.h
+++ b/src/lp_data/HighsLpSolverObject.h
@@ -22,13 +22,14 @@ class HighsLpSolverObject {
  public:
   HighsLpSolverObject(HighsLp& lp, HighsBasis& basis, HighsSolution& solution,
                       HighsInfo& highs_info, HEkk& ekk_instance,
-                      HighsCallback& callback, HighsOptions& options,
-                      HighsTimer& timer)
+                      HighsIpxStats& ipx_stats, HighsCallback& callback,
+                      HighsOptions& options, HighsTimer& timer)
       : lp_(lp),
         basis_(basis),
         solution_(solution),
         highs_info_(highs_info),
         ekk_instance_(ekk_instance),
+        ipx_stats_(ipx_stats),
         callback_(callback),
         options_(options),
         timer_(timer) {}
@@ -38,6 +39,7 @@ class HighsLpSolverObject {
   HighsSolution& solution_;
   HighsInfo& highs_info_;
   HEkk& ekk_instance_;
+  HighsIpxStats& ipx_stats_;
   HighsCallback& callback_;
   HighsOptions& options_;
   HighsTimer& timer_;
diff --git a/src/lp_data/HighsOptions.h b/src/lp_data/HighsOptions.h
index 82424a0adc..0543f547e5 100644
--- a/src/lp_data/HighsOptions.h
+++ b/src/lp_data/HighsOptions.h
@@ -333,6 +333,9 @@ struct HighsOptionsStruct {
 
   // Options for IPM solver
   HighsInt ipm_iteration_limit;
+  HighsInt cr1_iteration_limit;
+  HighsInt cr2_iteration_limit;
+  bool kkt_logging;
 
   // Options for PDLP solver
   bool pdlp_native_termination;
@@ -411,6 +414,7 @@ struct HighsOptionsStruct {
   // Options for MIP solver
   bool mip_detect_symmetry;
   bool mip_allow_restart;
+  bool mip_old_analytic_centre_method;
   HighsInt mip_max_nodes;
   HighsInt mip_max_stall_nodes;
   HighsInt mip_max_start_nodes;
@@ -421,6 +425,7 @@ struct HighsOptionsStruct {
   HighsInt mip_pool_soft_limit;
   HighsInt mip_pscost_minreliable;
   HighsInt mip_min_cliquetable_entries_for_parallelism;
+  HighsInt mip_compute_analytic_centre;
   HighsInt mip_report_level;
   double mip_feasibility_tolerance;
   double mip_rel_gap;
@@ -481,6 +486,9 @@ struct HighsOptionsStruct {
         output_flag(false),
         log_to_console(false),
         ipm_iteration_limit(0),
+        cr1_iteration_limit(0),
+        cr2_iteration_limit(0),
+        kkt_logging(false),
         pdlp_native_termination(false),
         pdlp_scaling(false),
         pdlp_iteration_limit(0),
@@ -545,6 +553,7 @@ struct HighsOptionsStruct {
         icrash_breakpoints(false),
         mip_detect_symmetry(false),
         mip_allow_restart(false),
+        mip_old_analytic_centre_method(false),
         mip_max_nodes(0),
         mip_max_stall_nodes(0),
         mip_max_start_nodes(0),
@@ -555,6 +564,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_report_level(0),
         mip_feasibility_tolerance(0.0),
         mip_rel_gap(0.0),
@@ -939,6 +949,12 @@ class HighsOptions : public HighsOptionsStruct {
                                        advanced, &mip_allow_restart, true);
     records.push_back(record_bool);
 
+    record_bool =
+        new OptionRecordBool("mip_old_analytic_centre_method",
+                             "Use origial or new analytic centre method in MIP",
+                             advanced, &mip_old_analytic_centre_method, true);
+    records.push_back(record_bool);
+
     record_int = new OptionRecordInt("mip_max_nodes",
                                      "MIP solver max number of nodes", advanced,
                                      &mip_max_nodes, 0, kHighsIInf, kHighsIInf);
@@ -1032,6 +1048,14 @@ class HighsOptions : public HighsOptionsStruct {
         kHighsIInf);
     records.push_back(record_int);
 
+    record_int = new OptionRecordInt(
+        "mip_compute_analytic_centre",
+        "Compute analytic centre for MIP: -1 => choose; 0 => off; 1 => on "
+        "(default)",
+        advanced, &mip_compute_analytic_centre, kMipAnalyticCentreCalulationMin,
+        kMipAnalyticCentreCalulationOn, kMipAnalyticCentreCalulationMax);
+    records.push_back(record_int);
+
     record_int =
         new OptionRecordInt("mip_report_level", "MIP solver reporting level",
                             now_advanced, &mip_report_level, 0, 1, 2);
@@ -1071,6 +1095,26 @@ class HighsOptions : public HighsOptionsStruct {
         &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf);
     records.push_back(record_int);
 
+    record_bool = new OptionRecordBool(
+        "kkt_logging",
+        "Perform logging in CR solver for KKT in IPX: Default = false",
+        advanced, &kkt_logging, false);
+    records.push_back(record_bool);
+
+    record_int =
+        new OptionRecordInt("cr1_iteration_limit",
+                            "Iteration limit for initial CR in IPX IPM solver: "
+                            "default is -1, so set internally",
+                            advanced, &cr1_iteration_limit, -1, -1, kHighsIInf);
+    records.push_back(record_int);
+
+    record_int =
+        new OptionRecordInt("cr2_iteration_limit",
+                            "Iteration limit for main CR in IPX IPM solver: "
+                            "default is -1, so set internally",
+                            advanced, &cr2_iteration_limit, -1, -1, kHighsIInf);
+    records.push_back(record_int);
+
     record_bool = new OptionRecordBool(
         "pdlp_native_termination",
         "Use native termination for PDLP solver: Default = false", advanced,
diff --git a/src/lp_data/HighsSolverStats.h b/src/lp_data/HighsSolverStats.h
new file mode 100644
index 0000000000..99b8d9f742
--- /dev/null
+++ b/src/lp_data/HighsSolverStats.h
@@ -0,0 +1,106 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/*                                                                       */
+/*    This file is part of the HiGHS linear optimization suite           */
+/*                                                                       */
+/*    Written and engineered 2008-2024 by Julian Hall, Ivet Galabova,    */
+/*    Leona Gottwald and Michael Feldmeier                               */
+/*                                                                       */
+/*    Available as open-source under the MIT License                     */
+/*                                                                       */
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/**@file lp_data/HighsSolverStats.h
+ * @brief Structs for HiGHS
+ */
+#ifndef LP_DATA_HIGHSSOLVERSTATS_H_
+#define LP_DATA_HIGHSSOLVERSTATS_H_
+
+#include <vector>
+
+#include "lp_data/HConst.h"
+
+enum HighsSolverStatsReport {
+  HighsSolverStatsReportPretty = 0,
+  HighsSolverStatsReportCsvHeader,
+  HighsSolverStatsReportCsvData
+};
+
+enum HighsSimplexWorkTerm {
+  HighsSimplexWorkTermInvertNumRow = 0,
+  HighsSimplexWorkTermInvertNumNz,
+  HighsSimplexWorkTermComputePD,
+  HighsSimplexWorkTermBtran,
+  HighsSimplexWorkTermPrice,
+  HighsSimplexWorkTermFtran,
+  HighsSimplexWorkTermFtranDse,
+  HighsSimplexWorkTermCount
+};
+
+const std::vector<std::string> kSimplexWorkNames = {
+    "InvertNumRow", "InvertNumNz", "ComputePD", "Btran",
+    "Price",        "Ftran",       "FtranDse"};
+const std::vector<double> kSimplexWorkCoefficients = {1.0, 2.0, 3.0, 4.0,
+                                                      5.0, 6.0, 7.0};
+
+enum HighsIpxWorkTerm {
+  HighsIpxWorkTermCr1IterNumRow = 0,
+  HighsIpxWorkTermCr1IterNumNz,
+  HighsIpxWorkTermCr2IterNumRow,
+  HighsIpxWorkTermCr2IterNumNz,
+  HighsIpxWorkTermCount
+};
+
+const std::vector<std::string> kIpxWorkNames = {
+    "CR1IterNumRow", "CR1IterNumNz", "CR2IterNumRow", "CR2IterNumNz"};
+const std::vector<double> kIpxWorkCoefficients = {1.0, 2.0, 3.0, 4.0};
+
+struct HighsSimplexStats {
+  bool valid;
+  HighsInt num_col;
+  HighsInt num_row;
+  HighsInt num_nz;
+  HighsInt iteration_count;
+  HighsInt num_invert;
+  HighsInt last_factored_basis_num_el;
+  HighsInt last_invert_num_el;
+  double col_aq_density;
+  double row_ep_density;
+  double row_ap_density;
+  double row_DSE_density;
+  double simplex_time;
+  void workTerms(double* terms) const;
+  double workEstimate() const;
+  void report(FILE* file, const std::string message = "",
+              const HighsInt style = HighsSolverStatsReportPretty) const;
+  void initialise(const HighsInt iteration_count_ = 0);
+};
+
+struct HighsIpxStats {
+  bool valid;
+  HighsInt num_col;
+  HighsInt num_row;
+  HighsInt num_nz;
+  HighsInt iteration_count;
+  std::vector<HighsInt> cr_type;
+  std::vector<HighsInt> cr_count;
+  std::vector<HighsInt> factored_basis_num_el;
+  std::vector<HighsInt> invert_num_el;
+  HighsInt num_type1_iteration;
+  HighsInt num_type2_iteration;
+  double average_type1_cr_count;
+  double average_type2_cr_count;
+  double average_type2_matrix_nz;
+  double average_type2_invert_nz;
+  double type1_time;
+  double basis0_time;
+  double type2_time;
+  double ipm_time;
+  double crossover_time;
+  void workTerms(double* terms);
+  double workEstimate();
+  void averages();
+  void report(FILE* file, const std::string message = "",
+              const HighsInt style = HighsSolverStatsReportPretty);
+  void initialise();
+};
+
+#endif /* LP_DATA_HIGHSSOLVERSTATS_H_ */
diff --git a/src/mip/HighsLpRelaxation.cpp b/src/mip/HighsLpRelaxation.cpp
index 4288f99437..73b675917c 100644
--- a/src/mip/HighsLpRelaxation.cpp
+++ b/src/mip/HighsLpRelaxation.cpp
@@ -1064,9 +1064,11 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
              int(lpsolver.getNumCol()), int(lpsolver.getNumRow()));
     }
   }
-  const bool solver_logging = false;
+  const bool solver_logging = false;  //! mipsolver.submip && !valid_basis;
   const bool detailed_simplex_logging = false;
-  if (solver_logging) lpsolver.setOptionValue("output_flag", true);
+  if (solver_logging) {
+    lpsolver.setOptionValue("output_flag", true);
+  }
   if (detailed_simplex_logging) {
     lpsolver.setOptionValue("output_flag", true);
     lpsolver.setOptionValue("log_dev_level", kHighsLogDevLevelVerbose);
@@ -1076,6 +1078,8 @@ HighsLpRelaxation::Status HighsLpRelaxation::run(bool resolve_on_error) {
 
   mipsolver.analysis_.mipTimerStart(simplex_solve_clock);
   HighsStatus callstatus = lpsolver.run();
+  lpsolver.setOptionValue("output_flag", false);  // !fix-2049
+  //  if (solver_logging) lpsolver.reportSimplexStats(stdout); // !fix-2049
   mipsolver.analysis_.mipTimerStop(simplex_solve_clock);
 
   const HighsInfo& info = lpsolver.getInfo();
diff --git a/src/mip/HighsMipSolver.cpp b/src/mip/HighsMipSolver.cpp
index 2b5b93e6e6..cf5ab3f7a2 100644
--- a/src/mip/HighsMipSolver.cpp
+++ b/src/mip/HighsMipSolver.cpp
@@ -134,7 +134,7 @@ void HighsMipSolver::run() {
   mipdata_->runPresolve(options_mip_->presolve_reduction_limit);
   analysis_.mipTimerStop(kMipClockRunPresolve);
   analysis_.mipTimerStop(kMipClockPresolve);
-  if (analysis_.analyse_mip_time & !submip)
+  if (analysis_.analyse_mip_time && !submip)
     highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                  "MIP-Timing: %11.2g - completed presolve\n", timer_.read());
   // Identify whether time limit has been reached (in presolve)
@@ -158,16 +158,53 @@ void HighsMipSolver::run() {
 
   analysis_.mipTimerStart(kMipClockSolve);
 
-  if (analysis_.analyse_mip_time & !submip)
+  if (analysis_.analyse_mip_time && !submip)
     highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                  "MIP-Timing: %11.2g - starting  setup\n", timer_.read());
   analysis_.mipTimerStart(kMipClockRunSetup);
   mipdata_->runSetup();
   analysis_.mipTimerStop(kMipClockRunSetup);
-  if (analysis_.analyse_mip_time & !submip)
+  if (analysis_.analyse_mip_time && !submip)
     highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                  "MIP-Timing: %11.2g - completed setup\n",
                  timer_.read(timer_.total_clock));
+  const bool initial_root_node_solve = false;
+  if (!submip && initial_root_node_solve) {
+    if (analysis_.analyse_mip_time) {
+      highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
+                   "MIP-Timing: %11.2g - starting  relaxation simplex solve\n",
+                   timer_.read());
+      analysis_.mipTimerStart(kMipClockRelaxationSimplexSolve);
+      analysis_.mipTimerStart(kMipClockSimplexNoBasisSolveLp);
+    }
+    Highs root;
+    HighsLp lpmodel(*model_);
+    lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0);
+    root.passModel(std::move(lpmodel));
+    root.setOptionValue("solve_relaxation", true);
+    root.setOptionValue("output_flag", true);
+    HighsStatus status = root.run();
+    HighsModelStatus root_modelstatus = root.getModelStatus();
+    if (status != HighsStatus::kOk) {
+      cleanupSolve();
+      return;
+    }
+    if (root_modelstatus != HighsModelStatus::kOptimal) {
+      modelstatus_ = root_modelstatus;
+      cleanupSolve();
+      return;
+    }
+    mipdata_->firstrootbasis = root.getBasis();
+    mipdata_->simplex_stats = root.getPresolvedLpSimplexStats();
+    mipdata_->simplex_stats.report(stdout, "Root node");
+    if (analysis_.analyse_mip_time) {
+      analysis_.mipTimerStop(kMipClockSimplexNoBasisSolveLp);
+      analysis_.mipTimerStop(kMipClockRelaxationSimplexSolve);
+      highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
+                   "MIP-Timing: %11.2g - completed relaxation simplex solve\n",
+                   timer_.read());
+    }
+  }
 restart:
   if (modelstatus_ == HighsModelStatus::kNotset) {
     // Check limits have not been reached before evaluating root node
@@ -186,7 +223,7 @@ void HighsMipSolver::run() {
       return;
     }
     analysis_.mipTimerStop(kMipClockTrivialHeuristics);
-    if (analysis_.analyse_mip_time & !submip)
+    if (analysis_.analyse_mip_time && !submip)
       highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                    "MIP-Timing: %11.2g - starting evaluate root node\n",
                    timer_.read(timer_.total_clock));
@@ -198,7 +235,7 @@ void HighsMipSolver::run() {
     if (analysis_.analyse_mip_time &&
         analysis_.mipTimerRunning(kMipClockIpmSolveLp))
       analysis_.mipTimerStop(kMipClockIpmSolveLp);
-    if (analysis_.analyse_mip_time & !submip)
+    if (analysis_.analyse_mip_time && !submip)
       highsLogUser(options_mip_->log_options, HighsLogType::kInfo,
                    "MIP-Timing: %11.2g - completed evaluate root node\n",
                    timer_.read(timer_.total_clock));
@@ -795,30 +832,43 @@ 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"
                "                    %.2f (solve)\n"
                "                    %.2f (postsolve)\n"
-               "  Max sub-MIP depth %d\n"
                "  Nodes             %llu\n"
                "  Repair LPs        %llu (%llu feasible; %llu iterations)\n"
                "  LP iterations     %llu (total)\n"
                "                    %llu (strong br.)\n"
                "                    %llu (separation)\n"
-               "                    %llu (heuristics)\n",
+               "                    %llu (heuristics)\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),
                analysis_.mipTimerRead(kMipClockPostsolve),
-               int(max_submip_level), (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);
+               (long long unsigned)mipdata_->heuristic_lp_iterations,
+               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/HighsMipSolver.h b/src/mip/HighsMipSolver.h
index 0a1eaae0ca..0a7a3395c6 100644
--- a/src/mip/HighsMipSolver.h
+++ b/src/mip/HighsMipSolver.h
@@ -92,10 +92,10 @@ class HighsMipSolver {
 
   ~HighsMipSolver();
 
-  void setModel(const HighsLp& model) {
-    model_ = &model;
-    solution_objective_ = kHighsInf;
-  }
+  //  void setModel(const HighsLp& model) {
+  //    model_ = &model;
+  //    solution_objective_ = kHighsInf;
+  //  }
 
   mutable HighsTimer timer_;
   void cleanupSolve();
diff --git a/src/mip/HighsMipSolverData.cpp b/src/mip/HighsMipSolverData.cpp
index 522804d18b..e3fee8fe8a 100644
--- a/src/mip/HighsMipSolverData.cpp
+++ b/src/mip/HighsMipSolverData.cpp
@@ -302,26 +302,88 @@ 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
+    const bool ac_logging = false;  //! mipsolver.submip;  // 2049 unset this
     Highs ipm;
-    ipm.setOptionValue("solver", "ipm");
-    ipm.setOptionValue("run_crossover", kHighsOffString);
+    if (mipsolver.options_mip_->mip_old_analytic_centre_method) {
+      // Original calculation is just IPM with crossover off
+      ipm.setOptionValue("solver", "ipm");
+      ipm.setOptionValue("run_crossover", kHighsOffString);
+    } else {
+      // True calculation performs centring properly
+      ipm.setOptionValue("run_centring", true);
+      ipm.setOptionValue("ipm_optimality_tolerance", 1e-2);
+    }
     ipm.setOptionValue("presolve", "off");
     ipm.setOptionValue("output_flag", false);
-    // ipm.setOptionValue("output_flag", !mipsolver.submip);
+    //    ipm.setOptionValue("output_flag", !mipsolver.submip);  // 2049 unset
+    //    this ultimately
     ipm.setOptionValue("ipm_iteration_limit", 200);
+    double time_available =
+        std::max(mipsolver.options_mip_->time_limit -
+                     mipsolver.timer_.read(mipsolver.timer_.total_clock),
+                 0.1);
+    ipm.setOptionValue("time_limit", time_available);
+    //    ipm.setOptionValue("kkt_logging", !mipsolver.submip);  // 2049 unset
+    //    this
+    // ultimately
+    //
+    // cr1_iteration_limit is what's set internal to IPX to limit the
+    // CR iterations before the initial basis is computed, and perhaps
+    // should not be changed, but maybe 500 is too high
+    HighsInt cr1_iteration_limit = 10 + mipsolver.model_->num_row_ / 20;
+    cr1_iteration_limit =
+        std::min(HighsInt(500), cr1_iteration_limit);  // IPX internal default
+    cr1_iteration_limit = std::min(
+        HighsInt(100), cr1_iteration_limit);  // 2049 set this ultimately
+                                              //  2049 set this ultimately
+    // ipm.setOptionValue("cr1_iteration_limit", cr1_iteration_limit);
+    // cr2_iteration_limit is not set internal to IPX, so the default
+    // value is m+100, which is OK if you're desperate to solve an LP,
+    // but can make the use of the AC in MIP prohibitively expensive
+    HighsInt cr2_iteration_limit = 50 + mipsolver.model_->num_row_ / 50;
+    cr2_iteration_limit = std::min(HighsInt(500), cr2_iteration_limit);
+    // 2049 Set this ultimately
+    // ipm.setOptionValue("cr2_iteration_limit", cr2_iteration_limit);
+
+    //    if (!mipsolver.submip) this->simplex_stats.report(stdout);2049 unset
+    //    this
+    // ultimately
+    //
+    ipm.passSimplexStats(this->simplex_stats);
     HighsLp lpmodel(*mipsolver.model_);
     lpmodel.col_cost_.assign(lpmodel.num_col_, 0.0);
     ipm.passModel(std::move(lpmodel));
 
-    //    if (!mipsolver.submip) {
-    //      const std::string file_name = mipsolver.model_->model_name_ +
-    //      "_ipm.mps"; printf("Calling ipm.writeModel(%s)\n",
-    //      file_name.c_str()); fflush(stdout); ipm.writeModel(file_name);
-    //    }
+    //    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);
+    double tt1 = mipsolver.timer_.read(mipsolver.timer_.total_clock);
+    if (ac_logging) {
+      HighsModelStatus ac_status = ipm.getModelStatus();
+      printf(
+          "grepAcLocal: model; num_row; CR limits+counts; time and status,"
+          "%s,%d,"
+          "%d,%d,"
+          "%d,%d,"
+          "%g,%s\n",
+          mipsolver.model_->model_name_.c_str(), int(ipm.getNumRow()),
+          int(cr1_iteration_limit), int(ipm.getInfo().max_cr_iteration_count1),
+          int(cr2_iteration_limit), int(ipm.getInfo().max_cr_iteration_count2),
+          tt1 - tt0, ipm.modelStatusToString(ac_status).c_str());
+      if (ac_status == HighsModelStatus::kSolveError ||
+          ac_status == HighsModelStatus::kUnknown) {
+        num_analytic_centre_fail++;
+      } else if (ac_status == HighsModelStatus::kOptimal) {
+        num_analytic_centre_opt++;
+      } else {
+        num_analytic_centre_other++;
+      }
+    }
     const std::vector<double>& sol = ipm.getSolution().col_value;
     if (HighsInt(sol.size()) != mipsolver.numCol()) return;
     analyticCenterStatus = ipm.getModelStatus();
@@ -345,6 +407,9 @@ void HighsMipSolverData::finishAnalyticCenterComputation(
     fflush(stdout);
   }
   analyticCenterComputed = true;
+  analyticCenterFailed =
+      analyticCenterStatus == HighsModelStatus::kSolveError ||
+      analyticCenterStatus == HighsModelStatus::kUnknown;
   if (analyticCenterStatus == HighsModelStatus::kOptimal) {
     HighsInt nfixed = 0;
     HighsInt nintfixed = 0;
@@ -381,6 +446,11 @@ void HighsMipSolverData::finishAnalyticCenterComputation(
                   nfixed, nintfixed);
     mipsolver.mipdata_->domain.propagate();
     if (mipsolver.mipdata_->domain.infeasible()) return;
+  } else if (!analyticCenterFailed) {
+    printf(
+        "HighsMipSolverData::finishAnalyticCenterComputation: "
+        "analyticCenterStatus = %s\n",
+        lp.getLpSolver().modelStatusToString(analyticCenterStatus).c_str());
   }
 }
 
@@ -611,6 +681,7 @@ void HighsMipSolverData::init() {
   firstlpsolobj = -kHighsInf;
   rootlpsolobj = -kHighsInf;
   analyticCenterComputed = false;
+  analyticCenterFailed = false;
   analyticCenterStatus = HighsModelStatus::kNotset;
   maxTreeSizeLog2 = 0;
   numRestarts = 0;
@@ -633,6 +704,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_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;
@@ -963,7 +1038,10 @@ void HighsMipSolverData::runSetup() {
   }
 #endif
 
-  if (upper_limit == kHighsInf) analyticCenterComputed = false;
+  if (upper_limit == kHighsInf) {
+    analyticCenterComputed = false;
+    analyticCenterFailed = false;
+  }
   analyticCenterStatus = HighsModelStatus::kNotset;
   analyticCenter.clear();
 
@@ -1851,8 +1929,6 @@ HighsLpRelaxation::Status HighsMipSolverData::evaluateRootLp() {
 }
 
 void HighsMipSolverData::evaluateRootNode() {
-  const bool compute_analytic_centre = true;
-  if (!compute_analytic_centre) printf("NOT COMPUTING ANALYTIC CENTRE!\n");
   HighsInt maxSepaRounds = mipsolver.submip ? 5 : kHighsIInf;
   if (numRestarts == 0)
     maxSepaRounds =
@@ -1865,10 +1941,15 @@ void HighsMipSolverData::evaluateRootNode() {
     startSymmetryDetection(tg, symData);
     mipsolver.analysis_.mipTimerStop(kMipClockStartSymmetryDetection);
   }
-  if (compute_analytic_centre && !analyticCenterComputed) {
-    highsLogUser(mipsolver.options_mip_->log_options, HighsLogType::kInfo,
-                 "MIP-Timing: %11.2g - starting analytic centre calculation\n",
-                 mipsolver.timer_.read());
+  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());
+      fflush(stdout);
+    }
     mipsolver.analysis_.mipTimerStart(kMipClockStartAnalyticCentreComputation);
     startAnalyticCenterComputation(tg);
     mipsolver.analysis_.mipTimerStop(kMipClockStartAnalyticCentreComputation);
@@ -2045,8 +2126,9 @@ void HighsMipSolverData::evaluateRootNode() {
       mipsolver.analysis_.mipTimerStop(kMipClockSeparation);
       return;
     }
-    if (nseparounds >= 5 && !mipsolver.submip && !analyticCenterComputed &&
-        compute_analytic_centre) {
+    if (nseparounds >= 5 && !mipsolver.submip &&
+        mipsolver.options_mip_->mip_compute_analytic_centre &&
+        !analyticCenterComputed) {
       if (checkLimits()) {
         mipsolver.analysis_.mipTimerStop(kMipClockSeparation);
         return;
@@ -2141,7 +2223,8 @@ void HighsMipSolverData::evaluateRootNode() {
   rootlpsolobj = lp.getObjective();
   lp.setIterationLimit(std::max(10000, int(10 * avgrootlpiters)));
 
-  if (!analyticCenterComputed && compute_analytic_centre) {
+  if (mipsolver.options_mip_->mip_compute_analytic_centre &&
+      !analyticCenterComputed) {
     if (checkLimits()) return;
 
     mipsolver.analysis_.mipTimerStart(kMipClockFinishAnalyticCentreComputation);
@@ -2263,7 +2346,8 @@ void HighsMipSolverData::evaluateRootNode() {
   if (lower_bound <= upper_limit) {
     if (!mipsolver.submip && mipsolver.options_mip_->mip_allow_restart &&
         mipsolver.options_mip_->presolve != kHighsOffString) {
-      if (!analyticCenterComputed && compute_analytic_centre) {
+      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 38693bf361..291bf6fba0 100644
--- a/src/mip/HighsMipSolverData.h
+++ b/src/mip/HighsMipSolverData.h
@@ -84,12 +84,14 @@ struct HighsMipSolverData {
   bool cliquesExtracted;
   bool rowMatrixSet;
   bool analyticCenterComputed;
+  bool analyticCenterFailed;
   HighsModelStatus analyticCenterStatus;
   bool detectSymmetries;
   HighsInt numRestarts;
   HighsInt numRestartsRoot;
   HighsInt numCliqueEntriesAfterPresolve;
   HighsInt numCliqueEntriesAfterFirstPresolve;
+  HighsSimplexStats simplex_stats;
 
   std::vector<HighsInt> ARstart_;
   std::vector<HighsInt> ARindex_;
@@ -119,6 +121,10 @@ struct HighsMipSolverData {
   HighsInt numintegercols;
   HighsInt maxTreeSizeLog2;
 
+  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;
@@ -170,6 +176,7 @@ struct HighsMipSolverData {
         cliquesExtracted(false),
         rowMatrixSet(false),
         analyticCenterComputed(false),
+        analyticCenterFailed(false),
         analyticCenterStatus(HighsModelStatus::kNotset),
         detectSymmetries(false),
         numRestarts(0),
@@ -184,6 +191,10 @@ struct HighsMipSolverData {
         rootlpsolobj(-kHighsInf),
         numintegercols(0),
         maxTreeSizeLog2(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 1c4aff5daf..851f01c700 100644
--- a/src/mip/HighsPrimalHeuristics.cpp
+++ b/src/mip/HighsPrimalHeuristics.cpp
@@ -132,8 +132,15 @@ bool HighsPrimalHeuristics::solveSubMip(
   submipoptions.presolve = "on";
   submipoptions.mip_detect_symmetry = false;
   submipoptions.mip_heuristic_effort = 0.8;
-  // setup solver and run it
+  // If the analytic centre calculation has failed for the parent MIP,
+  // then don't perform it for the sub-MIP
+
+  // 2049 Set this ultimately
+  //
+  if (mipsolver.mipdata_->analyticCenterFailed)
+    submipoptions.mip_compute_analytic_centre = kMipAnalyticCentreCalulationOff;
 
+  // setup solver and run it
   HighsSolution solution;
   solution.value_valid = false;
   solution.dual_valid = false;
diff --git a/src/mip/MipTimer.h b/src/mip/MipTimer.h
index bc6534cd93..fc927bdc42 100644
--- a/src/mip/MipTimer.h
+++ b/src/mip/MipTimer.h
@@ -24,6 +24,7 @@ enum iClockMip {
   kMipClockInit,
   kMipClockRunPresolve,
   kMipClockRunSetup,
+  kMipClockRelaxationSimplexSolve,
   kMipClockTrivialHeuristics,
   kMipClockEvaluateRootNode,
   kMipClockPerformAging0,
@@ -96,6 +97,8 @@ class MipTimer {
     clock[kMipClockInit] = timer_pointer->clock_def("Initialise");
     clock[kMipClockRunPresolve] = timer_pointer->clock_def("Run presolve");
     clock[kMipClockRunSetup] = timer_pointer->clock_def("Run setup");
+    clock[kMipClockRelaxationSimplexSolve] =
+        timer_pointer->clock_def("Relaxation simplex solve");
     clock[kMipClockTrivialHeuristics] =
         timer_pointer->clock_def("Trivial heuristics");
     clock[kMipClockEvaluateRootNode] =
@@ -255,6 +258,7 @@ class MipTimer {
     const std::vector<HighsInt> mip_clock_list{kMipClockInit,
                                                kMipClockRunPresolve,
                                                kMipClockRunSetup,
+                                               kMipClockRelaxationSimplexSolve,
                                                kMipClockTrivialHeuristics,
                                                kMipClockEvaluateRootNode,
                                                kMipClockPerformAging0,
diff --git a/src/simplex/HApp.h b/src/simplex/HApp.h
index b8de2a3f59..28cd8cf3f5 100644
--- a/src/simplex/HApp.h
+++ b/src/simplex/HApp.h
@@ -183,6 +183,7 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   // basis are taken from the unscaled solution of the scaled LP.
   bool solve_unscaled_lp = false;
   bool solved_unscaled_lp = false;
+  bool refine_solution = false;
   if (!incumbent_lp.scale_.has_scaling) {
     //
     // Solve the unscaled LP with unscaled NLA
@@ -206,7 +207,6 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
   } else {
     // Indicate that there is no (current) need to refine the solution
     // by solving the unscaled LP with scaled NLA
-    bool refine_solution = false;
     if (options.simplex_unscaled_solution_strategy ==
             kSimplexUnscaledSolutionStrategyNone ||
         options.simplex_unscaled_solution_strategy ==
@@ -334,22 +334,26 @@ inline HighsStatus solveLpSimplex(HighsLpSolverObject& solver_object) {
       if (ekk_instance.proofOfPrimalInfeasibility()) solve_unscaled_lp = false;
     }
     if (solve_unscaled_lp) {
+      if (refine_solution)
+        printf(
+            "GrepSolverStats: Model %s requires refinement due to num/max/sum "
+            "primal "
+            "(%" HIGHSINT_FORMAT "/%g/%g) and dual (%" HIGHSINT_FORMAT
+            "/%g/%g) "
+            "unscaled infeasibilities\n",
+            incumbent_lp.model_name_.c_str(),
+            highs_info.num_primal_infeasibilities,
+            highs_info.max_primal_infeasibility,
+            highs_info.sum_primal_infeasibilities,
+            highs_info.num_dual_infeasibilities,
+            highs_info.max_dual_infeasibility,
+            highs_info.sum_dual_infeasibilities);
       // Save options/strategies that may be changed
       HighsInt simplex_strategy = options.simplex_strategy;
       double dual_simplex_cost_perturbation_multiplier =
           options.dual_simplex_cost_perturbation_multiplier;
       HighsInt simplex_dual_edge_weight_strategy =
           ekk_info.dual_edge_weight_strategy;
-      // #1865 exposed that this should not be
-      // HighsModelStatus::kObjectiveBound, but
-      // HighsModelStatus::kObjectiveTarget, since if the latter is
-      // the model status for the scaled LP, any primal
-      // infeasibilities should be small, but must be cleaned up
-      // before (hopefully) a few phase 2 primal simplex iterations
-      // are required to attain the target for the unscaled LP
-      //
-      // In #1865, phase 2 primal simplex was forced with large primal
-      // infeasibilities
       if (num_unscaled_primal_infeasibilities == 0 ||
           scaled_model_status == HighsModelStatus::kObjectiveTarget) {
         // Only dual infeasibilities, or objective target reached (in
diff --git a/src/simplex/HEkk.cpp b/src/simplex/HEkk.cpp
index caa4e641e0..5006f2d275 100644
--- a/src/simplex/HEkk.cpp
+++ b/src/simplex/HEkk.cpp
@@ -3503,6 +3503,10 @@ HighsStatus HEkk::returnFromEkkSolve(const HighsStatus return_status) {
   // reverts to its value given by options_
   if (analysis_.analyse_simplex_time) analysis_.reportSimplexTimer();
   simplex_stats_.valid = true;
+  simplex_stats_.num_col = lp_.num_col_;
+  simplex_stats_.num_row = lp_.num_row_;
+  simplex_stats_.num_nz = lp_.a_matrix_.numNz();
+
   // Since HEkk::iteration_count_ includes iteration on presolved LP,
   // simplex_stats_.iteration_count is initialised to -
   // HEkk::iteration_count_
@@ -4421,29 +4425,97 @@ void HEkk::unitBtranResidual(const HighsInt row_out, const HVector& row_ep,
   }
 }
 
-void HighsSimplexStats::report(FILE* file, std::string message) const {
-  fprintf(file, "\nSimplex stats: %s\n", message.c_str());
-  fprintf(file, "   valid                      = %d\n", this->valid);
-  fprintf(file, "   iteration_count            = %d\n", this->iteration_count);
-  fprintf(file, "   num_invert                 = %d\n", this->num_invert);
-  fprintf(file, "   last_invert_num_el         = %d\n",
-          this->last_invert_num_el);
-  fprintf(file, "   last_factored_basis_num_el = %d\n",
-          this->last_factored_basis_num_el);
-  fprintf(file, "   col_aq_density             = %g\n", this->col_aq_density);
-  fprintf(file, "   row_ep_density             = %g\n", this->row_ep_density);
-  fprintf(file, "   row_ap_density             = %g\n", this->row_ap_density);
-  fprintf(file, "   row_DSE_density            = %g\n", this->row_DSE_density);
+void HEkk::passSimplexStats(const HighsSimplexStats simplex_stats) {
+  this->simplex_stats_ = simplex_stats;
+}
+
+void HighsSimplexStats::report(FILE* file, std::string message,
+                               const HighsInt style) const {
+  if (style == HighsSolverStatsReportPretty) {
+    fprintf(file, "\nSimplex stats: %s\n", message.c_str());
+    fprintf(file, "   valid                      = %d\n", this->valid);
+    fprintf(file, "   num_col                    = %d\n", this->num_col);
+    fprintf(file, "   num_row                    = %d\n", this->num_row);
+    fprintf(file, "   num_nz                     = %d\n", this->num_nz);
+    fprintf(file, "   iteration_count            = %d\n",
+            this->iteration_count);
+    fprintf(file, "   num_invert                 = %d\n", this->num_invert);
+    fprintf(file, "   last_factored_basis_num_el = %d\n",
+            this->last_factored_basis_num_el);
+    fprintf(file, "   last_invert_num_el         = %d\n",
+            this->last_invert_num_el);
+    fprintf(file, "   col_aq_density             = %g\n", this->col_aq_density);
+    fprintf(file, "   row_ep_density             = %g\n", this->row_ep_density);
+    fprintf(file, "   row_ap_density             = %g\n", this->row_ap_density);
+    fprintf(file, "   row_DSE_density            = %g\n",
+            this->row_DSE_density);
+    fprintf(file, "   simplex time =             = %g\n", this->simplex_time);
+  } else if (style == HighsSolverStatsReportCsvHeader) {
+    fprintf(file,
+            "valid,col,row,nz,iteration_count,num_invert,last_factored_basis_"
+            "num_el,last_invert_num_el,"
+            "col_aq_density,row_ep_density,row_ap_density,row_DSE_"
+            "density,simplex_time,");
+  } else if (style == HighsSolverStatsReportCsvData) {
+    fprintf(file, "%d,%d,%d,%d,%d,%d,%d,%d,%g,%g,%g,%g,%g,", int(this->valid),
+            int(this->num_col), int(this->num_row), int(this->num_nz),
+            int(this->iteration_count), int(this->num_invert),
+            int(this->last_factored_basis_num_el),
+            int(this->last_invert_num_el), this->col_aq_density,
+            this->row_ep_density, this->row_ap_density, this->row_DSE_density,
+            this->simplex_time);
+  } else {
+    fprintf(file, "Unknown simplex stats report style of %d\n", int(style));
+    assert(123 == 456);
+  }
 }
 
 void HighsSimplexStats::initialise(const HighsInt iteration_count_) {
   valid = false;
   iteration_count = -iteration_count_;
+  num_col = 0;
+  num_row = 0;
+  num_nz = 0;
   num_invert = 0;
-  last_invert_num_el = 0;
   last_factored_basis_num_el = 0;
+  last_invert_num_el = 0;
   col_aq_density = 0;
   row_ep_density = 0;
   row_ap_density = 0;
   row_DSE_density = 0;
+  simplex_time = 0;
+}
+
+void HighsSimplexStats::workTerms(double* terms) const {
+  const double nonbasic_nz =
+      double(this->num_nz + this->num_row - this->last_factored_basis_num_el);
+  terms[HighsSimplexWorkTermInvertNumRow] =
+      double(this->num_invert) * double(this->num_row);
+  terms[HighsSimplexWorkTermInvertNumNz] =
+      double(this->num_invert) * this->last_factored_basis_num_el;
+  terms[HighsSimplexWorkTermComputePD] =
+      double(this->num_invert) * double(this->last_invert_num_el + nonbasic_nz);
+  terms[HighsSimplexWorkTermBtran] = double(this->iteration_count) *
+                                     double(this->last_invert_num_el) *
+                                     this->row_ep_density;
+  terms[HighsSimplexWorkTermPrice] =
+      double(this->iteration_count) * nonbasic_nz * this->row_ep_density;
+  terms[HighsSimplexWorkTermFtran] = double(this->iteration_count) *
+                                     double(this->last_invert_num_el) *
+                                     this->col_aq_density;
+  terms[HighsSimplexWorkTermFtranDse] = double(this->iteration_count) *
+                                        double(this->last_invert_num_el) *
+                                        this->row_DSE_density;
+}
+
+double HighsSimplexStats::workEstimate() const {
+  double* terms = new double[HighsSimplexWorkTermCount];
+  this->workTerms(terms);
+  double work = 0;
+  for (HighsInt iX = 0; iX < HighsSimplexWorkTermCount; iX++) {
+    assert(terms[iX] > 0);
+    work += terms[iX] * kSimplexWorkCoefficients[iX];
+  }
+  delete[] terms;
+  return work;
 }
diff --git a/src/simplex/HEkk.h b/src/simplex/HEkk.h
index 0b253bebc0..e29226adc4 100644
--- a/src/simplex/HEkk.h
+++ b/src/simplex/HEkk.h
@@ -15,6 +15,7 @@
 #define SIMPLEX_HEKK_H_
 
 #include "lp_data/HighsCallback.h"
+#include "lp_data/HighsSolverStats.h"
 #include "simplex/HSimplexNla.h"
 #include "simplex/HighsSimplexAnalysis.h"
 #include "util/HSet.h"
@@ -157,9 +158,12 @@ class HEkk {
                            const vector<double>& rowUpper);
 
   const HighsSimplexStats& getSimplexStats() const { return simplex_stats_; }
+  void passSimplexStats(const HighsSimplexStats simplex_stats);
   void initialiseSimplexStats() { simplex_stats_.initialise(iteration_count_); }
-  void reportSimplexStats(FILE* file, const std::string message = "") const {
-    simplex_stats_.report(file, message);
+  void reportSimplexStats(
+      FILE* file, const std::string message = "",
+      const HighsInt style = HighsSolverStatsReportPretty) const {
+    simplex_stats_.report(file, message, style);
   }
 
   // Make this private later