diff --git a/cmake/sources-python.cmake b/cmake/sources-python.cmake
index ac1c94c7e4..e944f7b54c 100644
--- a/cmake/sources-python.cmake
+++ b/cmake/sources-python.cmake
@@ -387,6 +387,7 @@ set(highs_headers_python
     highs/presolve/ICrashUtil.h
     highs/presolve/ICrashX.h
     highs/presolve/PresolveComponent.h
+    highs/presolve/PresolveTimer.h
     highs/qpsolver/a_asm.hpp
     highs/qpsolver/a_quass.hpp
     highs/qpsolver/basis.hpp
diff --git a/cmake/sources.cmake b/cmake/sources.cmake
index 1e605de4b7..5b7ab6b305 100644
--- a/cmake/sources.cmake
+++ b/cmake/sources.cmake
@@ -544,6 +544,7 @@ set(highs_headers
     presolve/ICrashUtil.h
     presolve/ICrashX.h
     presolve/PresolveComponent.h
+    presolve/PresolveTimer.h
     qpsolver/a_asm.hpp
     qpsolver/a_quass.hpp
     qpsolver/basis.hpp
diff --git a/highs/io/HighsIO.cpp b/highs/io/HighsIO.cpp
index 7f0b45170a..648587c64b 100644
--- a/highs/io/HighsIO.cpp
+++ b/highs/io/HighsIO.cpp
@@ -314,6 +314,11 @@ const std::string highsBoolToString(const bool b, const HighsInt field_width) {
   return b ? " true" : "false";
 }
 
+const std::string highsIntToPlural(const HighsInt i, const bool y) {
+  if (y) return i == 1 ? "y" : "ies";
+  return i == 1 ? "" : "s";
+}
+
 const std::string highsTimeToString(const double time) {
   return
 #ifndef NDEBUG
diff --git a/highs/io/HighsIO.h b/highs/io/HighsIO.h
index c35ad89dd0..91a21b8275 100644
--- a/highs/io/HighsIO.h
+++ b/highs/io/HighsIO.h
@@ -111,6 +111,7 @@ std::string highsFormatToString(const char* format, ...);
 const std::string highsBoolToString(const bool b,
                                     const HighsInt field_width = 2);
 
+const std::string highsIntToPlural(const HighsInt i, const bool y = false);
 const std::string highsInsertMdEscapes(const std::string& from_string);
 const std::string highsInsertMdId(const std::string& from_string);
 const std::string highsTimeToString(const double time);
diff --git a/highs/lp_data/HConst.h b/highs/lp_data/HConst.h
index 32b9fc5007..5f8fff712f 100644
--- a/highs/lp_data/HConst.h
+++ b/highs/lp_data/HConst.h
@@ -92,12 +92,14 @@ enum HighsAnalysisLevel {
   kHighsAnalysisLevelNlaTime = 32,
   kHighsAnalysisLevelMipData = 64,
   kHighsAnalysisLevelMipTime = 128,
+  kHighsAnalysisLevelPresolveTime = 256,
   kHighsAnalysisLevelMin = kHighsAnalysisLevelNone,
   kHighsAnalysisLevelMax =
       kHighsAnalysisLevelModelData + kHighsAnalysisLevelSolverSummaryData +
       kHighsAnalysisLevelSolverRuntimeData + kHighsAnalysisLevelSolverTime +
       kHighsAnalysisLevelNlaData + kHighsAnalysisLevelNlaTime +
-      kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime
+      kHighsAnalysisLevelMipData + kHighsAnalysisLevelMipTime +
+      kHighsAnalysisLevelPresolveTime
 };
 
 enum class HighsVarType : uint8_t {
diff --git a/highs/mip/MipTimer.h b/highs/mip/MipTimer.h
index f1044de065..3d55feec73 100644
--- a/highs/mip/MipTimer.h
+++ b/highs/mip/MipTimer.h
@@ -11,7 +11,7 @@
 #ifndef MIP_MIPTIMER_H_
 #define MIP_MIPTIMER_H_
 
-// Clocks for profiling the MIP dual mip solver
+// Clocks for profiling the MIP solver
 enum iClockMip {
   kMipClockTotal = 0,
   kMipClockPresolve,
@@ -114,7 +114,7 @@ enum iClockMip {
   kNumMipClock  //!< Number of MIP clocks
 };
 
-const double tolerance_percent_report = 0.1;
+static const double kMipClockTolerancePercentReport = 0.1;
 
 class MipTimer {
  public:
@@ -371,7 +371,7 @@ class MipTimer {
                                                kMipClockSearch,
                                                kMipClockPostsolve};
     reportMipClockList("MipLevl1", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal, tolerance_percent_report);
+                       kMipClockTotal, kMipClockTolerancePercentReport);
   };
 
   void reportMipSolveLpClock(const HighsTimerClock& mip_timer_clock) {
@@ -385,20 +385,20 @@ class MipTimer {
         kMipClockHipoSolveLp,
         kMipClockIpxSolveLp};
     reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSubMipSolveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockSubMipSolve};
     reportMipClockList("MipSlvLp", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipPresolveClock(const HighsTimerClock& mip_timer_clock) {
     const std::vector<HighsInt> mip_clock_list{kMipClockProbingPresolve,
                                                kMipClockEnumerationPresolve};
     reportMipClockList("MipPrslv", mip_clock_list, mip_timer_clock,
-                       kMipClockRunPresolve, tolerance_percent_report);
+                       kMipClockRunPresolve, kMipClockTolerancePercentReport);
   };
 
   void reportAltEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) {
@@ -407,7 +407,7 @@ class MipTimer {
                                                kMipClockEvaluateRootNode2};
     reportMipClockList(
         "AltEvaluateRootNode", mip_clock_list, mip_timer_clock,
-        kMipClockEvaluateRootNode);  //, tolerance_percent_report);
+        kMipClockEvaluateRootNode);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipEvaluateRootNodeClock(const HighsTimerClock& mip_timer_clock) {
@@ -433,7 +433,7 @@ class MipTimer {
     };
     reportMipClockList(
         "MipEvaluateRootNode", mip_clock_list, mip_timer_clock,
-        kMipClockEvaluateRootNode);  //, tolerance_percent_report);
+        kMipClockEvaluateRootNode);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipRootSeparationClock(const HighsTimerClock& mip_timer_clock) {
@@ -442,8 +442,9 @@ class MipTimer {
         kMipClockRootSeparationFinishAnalyticCentreComputation,
         kMipClockRootSeparationCentralRounding,
         kMipClockRootSeparationEvaluateRootLp};
-    reportMipClockList("MipRootSeparation", mip_clock_list, mip_timer_clock,
-                       kMipClockRootSeparation);  //, tolerance_percent_report);
+    reportMipClockList(
+        "MipRootSeparation", mip_clock_list, mip_timer_clock,
+        kMipClockRootSeparation);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSearchClock(const HighsTimerClock& mip_timer_clock) {
@@ -455,7 +456,7 @@ class MipTimer {
         //	kMipClock@
     };
     reportMipClockList("MipSerch", mip_clock_list, mip_timer_clock,
-                       kMipClockSearch, tolerance_percent_report);
+                       kMipClockSearch, kMipClockTolerancePercentReport);
   };
 
   void reportMipDiveClock(const HighsTimerClock& mip_timer_clock) {
@@ -463,7 +464,7 @@ class MipTimer {
         kMipClockDiveEvaluateNode, kMipClockDivePrimalHeuristics,
         kMipClockTheDive, kMipClockBacktrackPlunge, kMipClockPerformAging2};
     reportMipClockList("MipDive_", mip_clock_list, mip_timer_clock,
-                       kMipClockDive, tolerance_percent_report);
+                       kMipClockDive, kMipClockTolerancePercentReport);
   };
 
   void reportMipDivePrimalHeuristicsClock(
@@ -472,7 +473,7 @@ class MipTimer {
         kMipClockDiveRandomizedRounding, kMipClockDiveRens, kMipClockDiveRins};
     reportMipClockList("MipDivePrimalHeuristics", mip_clock_list,
                        mip_timer_clock, kMipClockDivePrimalHeuristics,
-                       tolerance_percent_report);
+                       kMipClockTolerancePercentReport);
   };
 
   void reportMipNodeSearchClock(const HighsTimerClock& mip_timer_clock) {
@@ -481,8 +482,9 @@ class MipTimer {
         //      kMipClockSearchBacktrack,
         kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1,
         kMipClockNodeSearchSeparation};  //, kMipClockStoreBasis};
-    reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock,
-                       kMipClockNodeSearch);  //, tolerance_percent_report);
+    reportMipClockList(
+        "MipNodeSearch", mip_clock_list, mip_timer_clock,
+        kMipClockNodeSearch);  //, kMipClockTolerancePercentReport);
   };
 
   void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) {
@@ -490,7 +492,7 @@ class MipTimer {
         kMipClockImplboundSepa, kMipClockCliqueSepa, kMipClockTableauSepa,
         kMipClockPathAggrSepa, kMipClockModKSepa};
     reportMipClockList("MipSeparation", mip_clock_list, mip_timer_clock,
-                       kMipClockTotal);  //, tolerance_percent_report);
+                       kMipClockTotal);  //, kMipClockTolerancePercentReport);
   };
 
   void csvMipClock(const std::string model_name,
diff --git a/highs/pdlp/hipdlp/pdhg.cc b/highs/pdlp/hipdlp/pdhg.cc
index 1c603cf8ec..1f4cd604ae 100644
--- a/highs/pdlp/hipdlp/pdhg.cc
+++ b/highs/pdlp/hipdlp/pdhg.cc
@@ -1717,7 +1717,7 @@ void PDLPSolver::setup(const HighsOptions& options, HighsTimer& timer) {
 #ifdef CUPDLP_GPU
                "GPU"
 #else
-                "CPU: performance may be disappointing!"
+               "CPU: performance may be disappointing!"
 #endif
   );
 #ifdef CUPDLP_GPU
diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp
index 8f55b22be3..6885bebc43 100644
--- a/highs/presolve/HPresolve.cpp
+++ b/highs/presolve/HPresolve.cpp
@@ -25,6 +25,7 @@
 #include "mip/HighsObjectiveFunction.h"
 #include "mip/MipTimer.h"
 #include "presolve/HighsPostsolveStack.h"
+#include "presolve/PresolveTimer.h"
 #include "test_kkt/DevKkt.h"
 #include "util/HFactor.h"
 #include "util/HighsCDouble.h"
@@ -1739,12 +1740,13 @@ HPresolve::Result HPresolve::runProbing(HighsPostsolveStack& postsolve_stack) {
       // Check for timeout
       tt = this->timer->read();
       if (tt > options->time_limit) {
-        highsLogUser(options->log_options, HighsLogType::kInfo,
-                     "Time limit reached in probing: "
-                     "consider not using probing by setting option "
-                     "presolve_rule_off to 2^%-d = %d\n",
-                     int(kPresolveRuleProbing),
-                     int(std::pow(int(2), int(kPresolveRuleProbing))));
+        highsLogUser(
+            options->log_options, HighsLogType::kInfo,
+            "Time limit reached in probing: "
+            "consider not using probing by setting option "
+            "presolve_rule_off to 2^%-d = %d\n",
+            int(kPresolveRuleProbing),
+            int(std::pow(int(2), static_cast<int>(kPresolveRuleProbing))));
         return Result::kStopped;
       }
 
@@ -5664,13 +5666,16 @@ HPresolve::Result HPresolve::initialRowAndColPresolve(
   // arrays are not initialized, also unset changedRowFlag so that the row will
   // be added to the changed row vector when it is changed after it was
   // processed
+  analysis_.presolveTimerStart(kPresolveClockInitialRow);
   for (HighsInt row = 0; row != model->num_row_; ++row) {
     if (rowDeleted[row]) continue;
     HPRESOLVE_CHECKED_CALL(rowPresolve(postsolve_stack, row));
     changedRowFlag[row] = false;
   }
+  analysis_.presolveTimerStop(kPresolveClockInitialRow);
 
   // same for the columns
+  analysis_.presolveTimerStart(kPresolveClockInitialCol);
   for (HighsInt col = 0; col != model->num_col_; ++col) {
     if (colDeleted[col]) continue;
     // round and update bounds
@@ -5680,6 +5685,7 @@ HPresolve::Result HPresolve::initialRowAndColPresolve(
     HPRESOLVE_CHECKED_CALL(colPresolve(postsolve_stack, col));
     changedColFlag[col] = false;
   }
+  analysis_.presolveTimerStop(kPresolveClockInitialCol);
 
   return checkLimits(postsolve_stack);
 }
@@ -5689,15 +5695,25 @@ HPresolve::Result HPresolve::fastPresolveLoop(
   do {
     storeCurrentProblemSize();
 
+    analysis_.presolveTimerStart(kPresolveClockFastLoopRowSingletons);
     HPRESOLVE_CHECKED_CALL(removeRowSingletons(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockFastLoopRowSingletons);
 
+    analysis_.presolveTimerStart(kPresolveClockFastLoopColSingletons);
     HPRESOLVE_CHECKED_CALL(presolveChangedRows(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockFastLoopColSingletons);
 
+    analysis_.presolveTimerStart(kPresolveClockFastLoopDoubletonEquations);
     HPRESOLVE_CHECKED_CALL(removeDoubletonEquations(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockFastLoopDoubletonEquations);
 
+    analysis_.presolveTimerStart(kPresolveClockFastLoopChangedRows);
     HPRESOLVE_CHECKED_CALL(presolveColSingletons(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockFastLoopChangedRows);
 
+    analysis_.presolveTimerStart(kPresolveClockFastLoopChangedCols);
     HPRESOLVE_CHECKED_CALL(presolveChangedCols(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockFastLoopChangedCols);
 
   } while (problemSizeReduction() > 0.01);
 
@@ -5749,8 +5765,8 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
   // Set up the logic to allow presolve rules, and logging for their
   // effectiveness
   analysis_.setup(this->model, this->options, this->numDeletedRows,
-                  this->numDeletedCols, silent);
-
+                  this->numDeletedCols, silent, this->timer);
+  analysis_.presolveTimerStart(kPresolveClockPresolve);
   if (options->presolve != kHighsOffString) {
     if (mipsolver) mipsolver->mipdata_->cliquetable.setPresolveFlag(true);
 
@@ -5777,7 +5793,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
     assert(this->timer);
     assert(this->timer->running());
 
+    analysis_.presolveTimerStart(kPresolveClockInitial);
     HPRESOLVE_CHECKED_CALL(initialRowAndColPresolve(postsolve_stack));
+    analysis_.presolveTimerStop(kPresolveClockInitial);
 
     HighsInt numParallelRowColCalls = 0;
     // ReductionType::kEqualityRowAddition(s) has no basis postsolve,
@@ -5811,7 +5829,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
         report();
       }
 
+      analysis_.presolveTimerStart(kPresolveClockFastLoop);
       HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack));
+      analysis_.presolveTimerStop(kPresolveClockFastLoop);
 
       storeCurrentProblemSize();
 
@@ -5824,14 +5844,19 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
             applyConflictGraphSubstitutions(postsolve_stack, numDelCol));
       }
 
-      if (analysis_.allow_rule_[kPresolveRuleAggregator])
+      if (analysis_.allow_rule_[kPresolveRuleAggregator]) {
+        analysis_.presolveTimerStart(kPresolveClockAggregator);
         HPRESOLVE_CHECKED_CALL(aggregator(postsolve_stack));
+        analysis_.presolveTimerStop(kPresolveClockAggregator);
+      }
 
       if (problemSizeReduction() > 0.05) continue;
 
       if (trySparsify && analysis_.allow_rule_[kPresolveRuleSparsify]) {
         HighsInt numNz = numNonzeros();
+        analysis_.presolveTimerStart(kPresolveClockSparsify);
         HPRESOLVE_CHECKED_CALL(sparsify(postsolve_stack));
+        analysis_.presolveTimerStop(kPresolveClockSparsify);
         double nzReduction =
             100.0 * (1.0 - (numNonzeros() / static_cast<double>(numNz)));
 
@@ -5839,12 +5864,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
           highsLogDev(options->log_options, HighsLogType::kInfo,
                       "Sparsify removed %.1f%% of nonzeros\n", nzReduction);
 
-          // #1710 exposes that this should not be
-          //
-          // fastPresolveLoop(postsolve_stack);
-          //
-          // but
+          analysis_.presolveTimerStart(kPresolveClockFastLoop);
           HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack));
+          analysis_.presolveTimerStop(kPresolveClockFastLoop);
         }
         trySparsify = false;
       }
@@ -5853,7 +5875,11 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
           numParallelRowColCalls < 5) {
         if (shrinkProblemEnabled && (numDeletedCols >= model->num_col_ / 2 ||
                                      numDeletedRows >= model->num_row_ / 2)) {
+          //	analysis_.presolveTimerStart(kPresolveClock@);
+          //	analysis_.presolveTimerStop(kPresolveClock@);
+          analysis_.presolveTimerStart(kPresolveClockShrinkProblem);
           shrinkProblem(postsolve_stack);
+          analysis_.presolveTimerStop(kPresolveClockShrinkProblem);
 
           toCSC(model->a_matrix_.value_, model->a_matrix_.index_,
                 model->a_matrix_.start_);
@@ -5861,12 +5887,16 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
                     model->a_matrix_.start_);
         }
         storeCurrentProblemSize();
+        analysis_.presolveTimerStart(kPresolveClockParallelRowsAndCols);
         HPRESOLVE_CHECKED_CALL(detectParallelRowsAndCols(postsolve_stack));
+        analysis_.presolveTimerStop(kPresolveClockParallelRowsAndCols);
         ++numParallelRowColCalls;
         if (problemSizeReduction() > 0.05) continue;
       }
 
+      analysis_.presolveTimerStart(kPresolveClockFastLoop);
       HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack));
+      analysis_.presolveTimerStop(kPresolveClockFastLoop);
 
       if (mipsolver != nullptr) {
         HighsInt num_strengthened = -1;
@@ -5879,7 +5909,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
                       num_strengthened);
       }
 
+      analysis_.presolveTimerStart(kPresolveClockFastLoop);
       HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack));
+      analysis_.presolveTimerStop(kPresolveClockFastLoop);
 
       if (mipsolver != nullptr && numCliquesBeforeProbing == -1) {
         numCliquesBeforeProbing = mipsolver->mipdata_->cliquetable.numCliques();
@@ -5912,7 +5944,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
       if (!dependentEquationsCalled) {
         if (shrinkProblemEnabled && (numDeletedCols >= model->num_col_ / 2 ||
                                      numDeletedRows >= model->num_row_ / 2)) {
+          analysis_.presolveTimerStart(kPresolveClockShrinkProblem);
           shrinkProblem(postsolve_stack);
+          analysis_.presolveTimerStop(kPresolveClockShrinkProblem);
 
           toCSC(model->a_matrix_.value_, model->a_matrix_.index_,
                 model->a_matrix_.start_);
@@ -5921,11 +5955,15 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
         }
         storeCurrentProblemSize();
         if (analysis_.allow_rule_[kPresolveRuleDependentEquations]) {
+          analysis_.presolveTimerStart(kPresolveClockDependentEquations);
           HPRESOLVE_CHECKED_CALL(removeDependentEquations(postsolve_stack));
+          analysis_.presolveTimerStop(kPresolveClockDependentEquations);
           dependentEquationsCalled = true;
         }
         if (analysis_.allow_rule_[kPresolveRuleDependentFreeCols])
-          HPRESOLVE_CHECKED_CALL(removeDependentFreeCols(postsolve_stack));
+          analysis_.presolveTimerStart(kPresolveClockDependentFreeCol);
+        HPRESOLVE_CHECKED_CALL(removeDependentFreeCols(postsolve_stack));
+        analysis_.presolveTimerStop(kPresolveClockDependentFreeCol);
         if (problemSizeReduction() > 0.05) continue;
       }
 
@@ -5956,6 +5994,8 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) {
 
   if (mipsolver != nullptr) HPRESOLVE_CHECKED_CALL(scaleMIP(postsolve_stack));
 
+  analysis_.presolveTimerStop(kPresolveClockPresolve);
+  analysis_.reportPresolveTimer();
   // analysePresolveRuleLog() should return true - no errors
   assert(analysis_.analysePresolveRuleLog());
   // Possibly report presolve log
@@ -6397,19 +6437,30 @@ HPresolve::Result HPresolve::removeDependentEquations(
   const bool logging_on = analysis_.logging_on_;
   if (equations.empty()) return Result::kOk;
 
+  auto returnOk = [&]() {
+    analysis_.logging_on_ = logging_on;
+    if (logging_on)
+      analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations);
+    return Result::kOk;
+  };
+
   if (logging_on)
     analysis_.startPresolveRuleLog(kPresolveRuleDependentEquations);
 
   HighsSparseMatrix matrix;
-  matrix.num_col_ = equations.size();
+  HighsInt num_equations = equations.size();
+  matrix.num_col_ = num_equations;
   matrix.num_row_ = model->num_col_ + 1;
-  matrix.start_.resize(matrix.num_col_ + 1);
+  matrix.start_.resize(num_equations + 1);
   matrix.start_[0] = 0;
-  const HighsInt maxCapacity = numNonzeros() + matrix.num_col_;
+  const HighsInt maxCapacity = numNonzeros() + num_equations;
   matrix.value_.reserve(maxCapacity);
   matrix.index_.reserve(maxCapacity);
 
-  std::vector<HighsInt> eqSet(matrix.num_col_);
+  std::vector<HighsInt> eqSet(num_equations);
+  std::vector<HighsInt> row_count;
+  row_count.assign(model->num_col_, 0);
+
   HighsInt i = 0;
   for (const std::pair<HighsInt, HighsInt>& p : equations) {
     HighsInt eq = p.second;
@@ -6417,8 +6468,10 @@ HPresolve::Result HPresolve::removeDependentEquations(
 
     // add entries of equation
     for (const HighsSliceNonzero& nonz : getRowVector(eq)) {
+      HighsInt iCol = nonz.index();
+      row_count[iCol]++;
       matrix.value_.push_back(nonz.value());
-      matrix.index_.push_back(nonz.index());
+      matrix.index_.push_back(iCol);
     }
 
     // add entry for artificial rhs column
@@ -6429,7 +6482,22 @@ HPresolve::Result HPresolve::removeDependentEquations(
 
     matrix.start_[i] = matrix.value_.size();
   }
-  std::vector<HighsInt> colSet(matrix.num_col_);
+  // Find the number of (true) variables in the system of equations as
+  // the number of columns with entries in at least one equation
+  HighsInt num_variables = 0;
+  for (HighsInt iCol = 0; iCol < model->num_col_; iCol++)
+    if (row_count[iCol]) num_variables++;
+  HighsInt num_nz = matrix.numNz();
+  const bool silent = silentLog();
+  if (!silent)
+    highsLogUser(options->log_options, HighsLogType::kInfo,
+                 "Considering dependency of %d equation%s in %d variable%s "
+                 "with %d nonzero%s\n",
+                 int(num_equations), highsIntToPlural(num_equations).c_str(),
+                 int(num_variables), highsIntToPlural(num_variables).c_str(),
+                 int(num_nz), highsIntToPlural(num_nz).c_str());
+  // Identify any dependent equations
+  std::vector<HighsInt> colSet(num_equations);
   std::iota(colSet.begin(), colSet.end(), 0);
   HFactor factor;
   factor.setup(matrix, colSet);
@@ -6442,30 +6510,40 @@ HPresolve::Result HPresolve::removeDependentEquations(
   //
   // ToDo: This is strictly non-deterministic, but so conservative
   // that it'll only reap the cases when factor.build never finishes
-  const double time_limit =
-      std::max(1.0, std::min(0.01 * options->time_limit, 1000.0));
+  const double kMaxDependentEquationsTime = 100;
+  const double time_limit = std::max(
+      1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime));
   factor.setTimeLimit(time_limit);
-  const bool silent = silentLog();
   // Determine rank deficiency of the equations
   if (!silent)
     highsLogUser(options->log_options, HighsLogType::kInfo,
-                 "Dependent equations search running on %d equations with time "
+                 "Dependent equations search running with time "
                  "limit of %.2fs\n",
-                 static_cast<int>(matrix.num_col_), time_limit);
+                 time_limit);
   double time_taken = -this->timer->read();
   HighsInt build_return = factor.build();
   time_taken += this->timer->read();
+  // Analyse what's been removed
+  HighsInt num_removed_row = 0;
+  HighsInt num_removed_nz = 0;
+  HighsInt num_fictitious_rows_skipped = 0;
   if (build_return == kBuildKernelReturnTimeout) {
     // HFactor::build has timed out, so just return
-    if (!silent)
+    if (!silent) {
+      if (options->log_dev_level > 0)
+        highsLogUser(
+            options->log_options, HighsLogType::kInfo,
+            "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g,Terminated\n",
+            model->model_name_.c_str(), static_cast<int>(num_equations),
+            static_cast<int>(num_variables), static_cast<int>(model->num_col_),
+            static_cast<int>(num_nz), static_cast<int>(num_removed_row),
+            static_cast<int>(num_removed_nz), time_taken);
       highsLogUser(options->log_options, HighsLogType::kInfo,
                    "Dependent equations search terminated after %.3gs due to "
                    "expected time exceeding limit\n",
                    time_taken);
-    analysis_.logging_on_ = logging_on;
-    if (logging_on)
-      analysis_.stopPresolveRuleLog(kPresolveRuleDependentFreeCols);
-    return Result::kOk;
+    }
+    return returnOk();
   } else {
     double pct_off_timeout =
         1e2 * std::fabs(time_taken - time_limit) / time_limit;
@@ -6479,10 +6557,6 @@ HPresolve::Result HPresolve::removeDependentEquations(
   // build_return as rank_deficiency must be valid
   assert(build_return >= 0);
   const HighsInt rank_deficiency = build_return;
-  // Analyse what's been removed
-  HighsInt num_removed_row = 0;
-  HighsInt num_removed_nz = 0;
-  HighsInt num_fictitious_rows_skipped = 0;
   for (HighsInt k = 0; k < rank_deficiency; k++) {
     if (factor.var_with_no_pivot[k] >= 0) {
       HighsInt redundant_row = eqSet[factor.var_with_no_pivot[k]];
@@ -6494,22 +6568,29 @@ HPresolve::Result HPresolve::removeDependentEquations(
       num_fictitious_rows_skipped++;
     }
   }
-  if (!silent)
-    highsLogUser(options->log_options, HighsLogType::kInfo,
-                 "Dependent equations search removed %d rows and %d nonzeros "
-                 "in %.2fs (limit = %.2fs)\n",
-                 static_cast<int>(num_removed_row),
-                 static_cast<int>(num_removed_nz), time_taken, time_limit);
-  if (num_fictitious_rows_skipped)
-    highsLogDev(options->log_options, HighsLogType::kInfo,
-                ", avoiding %d fictitious rows",
-                static_cast<int>(num_fictitious_rows_skipped));
-  highsLogDev(options->log_options, HighsLogType::kInfo, "\n");
-
-  analysis_.logging_on_ = logging_on;
-  if (logging_on)
-    analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations);
-  return Result::kOk;
+  if (!silent) {
+    highsLogUser(
+        options->log_options, HighsLogType::kInfo,
+        "Search of %d equation%s with %d / %d variable%s and %d nonzero%s "
+        "removed %d dependent equation%s and %d nonzero%s "
+        "in %.2fs with bounds in (%.2fs, %.2fs) and limit = %.2fs)",
+        // clang-format off
+        static_cast<int>(num_equations), highsIntToPlural(num_equations).c_str(),
+        static_cast<int>(num_variables),
+	static_cast<int>(model->num_col_), highsIntToPlural(num_variables).c_str(),
+	static_cast<int>(num_nz), highsIntToPlural(num_nz).c_str(),
+	static_cast<int>(num_removed_row), highsIntToPlural(num_removed_row).c_str(),
+        static_cast<int>(num_removed_nz), highsIntToPlural(num_removed_nz).c_str(),
+        // clang-format on
+        time_taken, factor.min_time_bound_, factor.max_time_bound_, time_limit);
+    if (num_fictitious_rows_skipped)
+      highsLogDev(options->log_options, HighsLogType::kInfo,
+                  ", avoiding %d fictitious row%s",
+                  static_cast<int>(num_fictitious_rows_skipped),
+                  highsIntToPlural(num_fictitious_rows_skipped).c_str());
+    highsLogUser(options->log_options, HighsLogType::kInfo, "\n");
+  }
+  return returnOk();
 }
 
 HPresolve::Result HPresolve::removeDependentFreeCols(
diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp
index f4724ebd05..8ca3362847 100644
--- a/highs/presolve/HPresolveAnalysis.cpp
+++ b/highs/presolve/HPresolveAnalysis.cpp
@@ -7,34 +7,57 @@
 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 #include "lp_data/HighsModelUtils.h"
 #include "presolve/HPresolve.h"
+#include "presolve/PresolveTimer.h"
 
 void HPresolveAnalysis::setup(const HighsLp* model_,
                               const HighsOptions* options_,
                               const HighsInt& numDeletedRows_,
                               const HighsInt& numDeletedCols_,
-                              const bool silent) {
+                              const bool silent, HighsTimer* timer) {
   model = model_;
   options = options_;
   numDeletedRows = &numDeletedRows_;
   numDeletedCols = &numDeletedCols_;
 
+  timer_ = timer;
+  analyse_presolve_time_ =
+      kHighsAnalysisLevelPresolveTime & options->highs_analysis_level &&
+      !model_->isMip();
+  if (analyse_presolve_time_) {
+    HighsTimerClock clock;
+    clock.timer_pointer_ = timer_;
+    PresolveTimer presolve_timer;
+    presolve_timer.initialisePresolveClocks(clock);
+    presolve_clocks_ = clock;
+  }
+
   this->allow_rule_.assign(kPresolveRuleCount, true);
 
-  if (options->presolve_rule_off || options_->log_dev_level) {
+  if (!silent && options_->log_dev_level) {
+    // State which rules can be off, and what bit to set
+    highsLogUser(options->log_options, HighsLogType::kInfo,
+                 "Permitted suppression of presolve rules via "
+                 "presolve_rule_off option:\n");
+    HighsInt bit =
+        std::pow(int(2), static_cast<int>(kPresolveRuleFirstAllowOff));
+    for (HighsInt rule_type = kPresolveRuleFirstAllowOff;
+         rule_type < kPresolveRuleCount; rule_type++) {
+      // This is a rule that can be switched off
+      highsLogUser(options->log_options, HighsLogType::kInfo,
+                   "   Rule %2d (set bit %2d = %5d): %s\n", int(rule_type),
+                   int(rule_type), int(bit),
+                   utilPresolveRuleTypeToString(rule_type).c_str());
+      bit *= 2;
+    }
+  }
+  if (options->presolve_rule_off) {
     // Some presolve rules are off
     //
     // Transform options->presolve_rule_off into logical settings in
     // allow_rule_[*], commenting on the rules switched off
-    if (!silent) {
-      if (options->presolve_rule_off) {
-        highsLogUser(options->log_options, HighsLogType::kInfo,
-                     "Presolve rules not allowed:\n");
-      } else {
-        highsLogUser(options->log_options, HighsLogType::kInfo,
-                     "Permitted suppression of presolve rules via "
-                     "presolve_rule_off option:\n");
-      }
-    }
+    if (!silent)
+      highsLogUser(options->log_options, HighsLogType::kInfo,
+                   "Presolve rules not allowed:\n");
     HighsInt bit = 1;
     for (HighsInt rule_type = kPresolveRuleMin; rule_type < kPresolveRuleCount;
          rule_type++) {
@@ -44,13 +67,11 @@ void HPresolveAnalysis::setup(const HighsLp* model_,
         // This is a rule that can be switched off, so comment
         // positively if it is off
         allow_rule_[rule_type] = allow;
-        if (!silent)
-          if (!allow ||
-              (!options->presolve_rule_off && options_->log_dev_level))
-            highsLogUser(options->log_options, HighsLogType::kInfo,
-                         "   Rule %2d (set bit %2d = %5d): %s\n",
-                         int(rule_type), int(rule_type), int(bit),
-                         utilPresolveRuleTypeToString(rule_type).c_str());
+        if (!allow && !silent)
+          highsLogUser(options->log_options, HighsLogType::kInfo,
+                       "   Rule %2d (set bit %2d = %5d): %s\n", int(rule_type),
+                       int(rule_type), int(bit),
+                       utilPresolveRuleTypeToString(rule_type).c_str());
       } else if (!allow && !silent) {
         // This is a rule that cannot be switched off so, if an
         // attempt is made, don't allow it to be off and comment
@@ -237,3 +258,22 @@ bool HPresolveAnalysis::analysePresolveRuleLog(const bool report) {
   }
   return true;
 }
+
+void HPresolveAnalysis::presolveTimerStart(
+    const HighsInt presolve_clock) const {
+  if (!analyse_presolve_time_) return;
+  HighsInt highs_timer_clock = presolve_clocks_.clock_[presolve_clock];
+  presolve_clocks_.timer_pointer_->start(highs_timer_clock);
+}
+
+void HPresolveAnalysis::presolveTimerStop(const HighsInt presolve_clock) const {
+  if (!analyse_presolve_time_) return;
+  HighsInt highs_timer_clock = presolve_clocks_.clock_[presolve_clock];
+  presolve_clocks_.timer_pointer_->stop(highs_timer_clock);
+}
+
+void HPresolveAnalysis::reportPresolveTimer() {
+  if (!analyse_presolve_time_) return;
+  PresolveTimer presolve_timer;
+  presolve_timer.reportPresolveCoreClock(model->model_name_, presolve_clocks_);
+}
diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h
index 9fa48f32d2..e61eec96ca 100644
--- a/highs/presolve/HPresolveAnalysis.h
+++ b/highs/presolve/HPresolveAnalysis.h
@@ -8,10 +8,16 @@
 /**@file presolve/HPresolveAnalysis.h
  * @brief
  */
-#ifndef PRESOLVE_HIGHS_PRESOLVE_ANALYSIS_H_
-#define PRESOLVE_HIGHS_PRESOLVE_ANALYSIS_H_
+#ifndef PRESOLVE_HPRESOLVEANALYSIS_H_
+#define PRESOLVE_HPRESOLVEANALYSIS_H_
+
+#include "util/HighsTimer.h"
 
 class HPresolveAnalysis {
+ public:
+  HPresolveAnalysis() : timer_(nullptr), analyse_presolve_time_(false) {}
+
+  HighsTimer* timer_;
   const HighsLp* model;
   const HighsOptions* options;
   const bool* allow_rule;
@@ -33,20 +39,27 @@ class HPresolveAnalysis {
   HighsInt num_deleted_cols0_;
   HighsPresolveLog presolve_log_;
 
+  HighsTimerClock presolve_clocks_;
+  bool analyse_presolve_time_;
+
   // for LP presolve
   //
   // Transform options->presolve_rule_off into logical settings in
   // allow_rule_[*], commenting on the rules switched off
   void setup(const HighsLp* model_, const HighsOptions* options_,
              const HighsInt& numDeletedRows_, const HighsInt& numDeletedCols_,
-             const bool silent);
+             const bool silent, HighsTimer* timer);
+  void setupPresolveTime(const HighsOptions& options);
   void resetNumDeleted();
 
   std::string presolveReductionTypeToString(const HighsInt reduction_type);
   void startPresolveRuleLog(const HighsInt rule_type);
   void stopPresolveRuleLog(const HighsInt rule_type);
   bool analysePresolveRuleLog(const bool report = false);
+  void presolveTimerStart(const HighsInt presolve_clock = 0) const;
+  void presolveTimerStop(const HighsInt presolve_clock = 0) const;
+  void reportPresolveTimer();
   friend class HPresolve;
 };
 
-#endif
+#endif /* PRESOLVE_HPRESOLVEANALYSIS_H_ */
diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h
new file mode 100644
index 0000000000..138c6fb418
--- /dev/null
+++ b/highs/presolve/PresolveTimer.h
@@ -0,0 +1,159 @@
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/*                                                                       */
+/*    This file is part of the HiGHS linear optimization suite           */
+/*                                                                       */
+/*    Available as open-source under the MIT License                     */
+/*                                                                       */
+/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
+/**@file presolve/PresolveTimer.h
+ * @brief Indices of presolve iClocks
+ */
+#ifndef PRESOLVE_PRESOLVETIMER_H_
+#define PRESOLVE_PRESOLVETIMER_H_
+
+// Clocks for profiling presolve
+enum iClockPresolve {
+  kPresolveClockTotal = 0,
+  kPresolveClockPresolve,
+  kPresolveClockInitial,
+  kPresolveClockInitialRow,
+  kPresolveClockInitialCol,
+  kPresolveClockFastLoop,
+  kPresolveClockFastLoopRowSingletons,
+  kPresolveClockFastLoopColSingletons,
+  kPresolveClockFastLoopDoubletonEquations,
+  kPresolveClockFastLoopChangedRows,
+  kPresolveClockFastLoopChangedCols,
+  kPresolveClockAggregator,
+  kPresolveClockSparsify,
+  kPresolveClockParallelRowsAndCols,
+  kPresolveClockDependentEquations,
+  kPresolveClockDependentFreeCol,
+  kPresolveClockShrinkProblem,
+  //  kPresolveClock@,
+  kNumPresolveClock  //!< Number of PRESOLVE clocks
+};
+
+static const double kPresolveClockTolerancePercentReport = 0.1;
+
+class PresolveTimer {
+ public:
+  void initialisePresolveClocks(HighsTimerClock& presolve_timer_clock) {
+    HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_;
+    std::vector<HighsInt>& clock = presolve_timer_clock.clock_;
+
+    clock.resize(kNumPresolveClock);
+    clock[kPresolveClockTotal] = 0;
+    clock[kPresolveClockPresolve] = timer_pointer->clock_def("Presolve");
+    clock[kPresolveClockInitial] = timer_pointer->clock_def("Initial");
+    clock[kPresolveClockInitialRow] = timer_pointer->clock_def("Initial row");
+    clock[kPresolveClockInitialCol] = timer_pointer->clock_def("Initial col");
+    clock[kPresolveClockFastLoop] = timer_pointer->clock_def("Fast loop");
+    clock[kPresolveClockFastLoopRowSingletons] =
+        timer_pointer->clock_def("Fast loop: row singletons");
+    clock[kPresolveClockFastLoopColSingletons] =
+        timer_pointer->clock_def("Fast loop: col singletons");
+    clock[kPresolveClockFastLoopDoubletonEquations] =
+        timer_pointer->clock_def("Fast loop: doubleton equations");
+    clock[kPresolveClockFastLoopChangedRows] =
+        timer_pointer->clock_def("Fast loop: changed rows");
+    clock[kPresolveClockFastLoopChangedCols] =
+        timer_pointer->clock_def("Fast loop: changed cols");
+    clock[kPresolveClockAggregator] = timer_pointer->clock_def("Aggregator");
+    clock[kPresolveClockSparsify] = timer_pointer->clock_def("Sparsify");
+    clock[kPresolveClockParallelRowsAndCols] =
+        timer_pointer->clock_def("Parallel rows and cols");
+    clock[kPresolveClockDependentEquations] =
+        timer_pointer->clock_def("Dependent equations");
+    clock[kPresolveClockDependentFreeCol] =
+        timer_pointer->clock_def("Dependent free columns");
+    clock[kPresolveClockShrinkProblem] =
+        timer_pointer->clock_def("Shrink problem");
+    //    clock[kPresolveClock@] = timer_pointer->clock_def("@");
+  };
+
+  bool reportPresolveClockList(
+      const char* grepStamp, const std::vector<HighsInt> presolve_clock_list,
+      const HighsTimerClock& presolve_timer_clock,
+      const HighsInt kPresolveClockIdeal = kPresolveClockPresolve,
+      const double tolerance_percent_report_ = -1) {
+    HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_;
+    if (!timer_pointer->printf_flag) return false;
+    const std::vector<HighsInt>& clock = presolve_timer_clock.clock_;
+    HighsInt presolve_clock_list_size = presolve_clock_list.size();
+    std::vector<HighsInt> clockList;
+    clockList.resize(presolve_clock_list_size);
+    for (HighsInt en = 0; en < presolve_clock_list_size; en++) {
+      clockList[en] = clock[presolve_clock_list[en]];
+    }
+    const double ideal_sum_time =
+        timer_pointer->clock_time[clock[kPresolveClockIdeal]];
+    const double tolerance_percent_report =
+        tolerance_percent_report_ >= 0 ? tolerance_percent_report_ : 1e-8;
+    return timer_pointer->reportOnTolerance(
+        grepStamp, clockList, ideal_sum_time, tolerance_percent_report);
+  };
+
+  void csvPresolveClockList(const std::string& grep_query,
+                            const std::string& model_name,
+                            const std::vector<HighsInt> presolve_clock_list,
+                            const HighsTimerClock& presolve_timer_clock,
+                            const HighsInt kPresolveClockIdeal,
+                            const bool header, const bool end_line) {
+    HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_;
+    if (!timer_pointer->printf_flag) return;
+    const std::vector<HighsInt>& clock = presolve_timer_clock.clock_;
+    const double ideal_sum_time =
+        timer_pointer->clock_time[clock[kPresolveClockIdeal]];
+    if (ideal_sum_time < 1e-2) return;
+    const HighsInt num_clock = presolve_clock_list.size();
+    if (header) {
+      printf("grep_%s,model,ideal", grep_query.c_str());
+      for (HighsInt iX = 0; iX < num_clock; iX++) {
+        HighsInt iclock = clock[presolve_clock_list[iX]];
+        printf(",%s", timer_pointer->clock_names[iclock].c_str());
+      }
+      printf(",Unaccounted");
+      if (end_line) printf("\n");
+      return;
+    }
+    double sum_time = 0;
+    printf("grep_%s,%s,%11.4g", grep_query.c_str(), model_name.c_str(),
+           ideal_sum_time);
+    for (HighsInt iX = 0; iX < num_clock; iX++) {
+      HighsInt iclock = clock[presolve_clock_list[iX]];
+      double time = timer_pointer->read(iclock);
+      sum_time += time;
+      printf(",%11.4g", time);
+    }
+    printf(",%11.4g", ideal_sum_time - sum_time);
+    if (end_line) printf("\n");
+  }
+
+  void reportPresolveCoreClock(const std::string& model_name,
+                               const HighsTimerClock& presolve_timer_clock) {
+    const std::vector<HighsInt> presolve_clock_list{
+        //      kPresolveClockInitial,
+        kPresolveClockInitialRow, kPresolveClockInitialCol,
+        //	kPresolveClockFastLoop,
+        kPresolveClockFastLoopRowSingletons,
+        kPresolveClockFastLoopColSingletons,
+        kPresolveClockFastLoopDoubletonEquations,
+        kPresolveClockFastLoopChangedRows, kPresolveClockFastLoopChangedCols,
+        kPresolveClockAggregator, kPresolveClockSparsify,
+        kPresolveClockParallelRowsAndCols, kPresolveClockDependentEquations,
+        kPresolveClockDependentFreeCol, kPresolveClockShrinkProblem
+        //	kPresolveClock@
+    };
+    reportPresolveClockList("PresolveCore_", presolve_clock_list,
+                            presolve_timer_clock, kPresolveClockPresolve, 0.1);
+    csvPresolveClockList("GrepPresolveCore_", model_name, presolve_clock_list,
+                         presolve_timer_clock, kPresolveClockPresolve, true,
+                         true);
+    csvPresolveClockList("GrepPresolveCore_", model_name, presolve_clock_list,
+                         presolve_timer_clock, kPresolveClockPresolve, false,
+                         true);
+  };
+};
+
+#endif /* PRESOLVE_PRESOLVETIMER_H_ */
diff --git a/highs/simplex/HighsSimplexAnalysis.cpp b/highs/simplex/HighsSimplexAnalysis.cpp
index 1d24b7b2fb..1e9f1aba80 100644
--- a/highs/simplex/HighsSimplexAnalysis.cpp
+++ b/highs/simplex/HighsSimplexAnalysis.cpp
@@ -1227,10 +1227,10 @@ void HighsSimplexAnalysis::updateInvertFormData(const HFactor& factor) {
 
     HighsInt kernel_invert_num_el =
         factor.invert_num_el -
-        (factor.basis_matrix_num_el - factor.kernel_num_el);
+        (factor.basis_matrix_num_el - factor.kernel_num_el) - factor.kernel_dim;
     assert(factor.kernel_num_el);
-    double kernel_fill_factor =
-        (1.0 * kernel_invert_num_el) / factor.kernel_num_el;
+    double kernel_fill_factor = (1.0 * kernel_invert_num_el) /
+                                (factor.kernel_num_el + factor.kernel_dim);
     sum_kernel_fill_factor += kernel_fill_factor;
     running_average_kernel_fill_factor =
         0.95 * running_average_kernel_fill_factor + 0.05 * kernel_fill_factor;
diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp
index 31a16af254..38b02ccc8d 100644
--- a/highs/util/HFactor.cpp
+++ b/highs/util/HFactor.cpp
@@ -677,7 +677,11 @@ void HFactor::buildSimple() {
     b_start[iCol + 1] = BcountX;
     b_var[iCol] = iMat;
   }
-  // Record the number of elements in the basis matrix
+  // Record the number of elements in the basis matrix, remebering
+  // that BcountX is the number of entries in the nwork structural
+  // columns, so have to add num_row - nwork to get the entries in
+  // logical columns. In particular, if the basis matrix is an
+  // identity, BcountX = nwork = 0
   basis_matrix_num_el = num_row - nwork + BcountX;
 
   // count1 = 0;
@@ -819,16 +823,20 @@ void HFactor::buildSimple() {
   row_link_first.assign(num_basic + 1, -1);
   mr_count.assign(num_row, 0);
   HighsInt mr_countX = 0;
-  // Determine the number of entries in the kernel
-  kernel_num_el = 0;
+  // Determine the initial number of active nonzeros - to be updated
+  // as values are eliminated and fill-in or cancellation occur
+  num_active_nz_ = 0;
+  HighsInt check_nwork = 0;
   for (HighsInt iRow = 0; iRow < num_row; iRow++) {
     HighsInt count = mr_count_before[iRow];
     if (count > 0) {
+      // In the active part of the kernel
       mr_start[iRow] = mr_countX;
       mr_space[iRow] = count * 2;
       mr_countX += count * 2;
       rlinkAdd(iRow, count);
-      kernel_num_el += count + 1;
+      num_active_nz_ += count;
+      check_nwork++;
     }
   }
   mr_index.resize(mr_countX);
@@ -853,9 +861,11 @@ void HFactor::buildSimple() {
       const HighsInt iRow = b_index[k];
       const double value = b_value[k];
       if (mr_count_before[iRow] > 0) {
+        // In the active part of the kernel
         colInsert(iCol, iRow, value);
         rowInsert(iCol, iRow);
       } else {
+        // Above the active part of the kernel
         colStoreN(iCol, iRow, value);
       }
     }
@@ -863,8 +873,11 @@ void HFactor::buildSimple() {
     clinkAdd(iCol, mc_count_a[iCol]);
   }
   build_synthetic_tick += (num_row + nwork + MCcountX) * 40 + mr_countX * 20;
-  // Record the kernel dimension
+  // Record the dimension and number of entries in the kernel.
   kernel_dim = nwork;
+  kernel_num_el = num_active_nz_;
+  min_time_bound_ = kHighsInf;
+  max_time_bound_ = 0;
   assert((HighsInt)this->refactor_info_.pivot_row.size() == num_basic - nwork);
 }
 
@@ -875,40 +888,124 @@ HighsInt HFactor::buildKernel() {
   double fake_fill = 0;
   double fake_eliminate = 0;
 
+  HighsInt search_k = 0;
   const bool progress_report = false;  // num_basic != num_row;
   const HighsInt progress_frequency = 10000;
+
+  // Normally this->time_limit_ is kHighsInf, but if there is a finite
+  // time limit, it implies that buildKernel can bail out. This is
+  // (currently) used only with the dependent equations rule in
+  // presolve
+  //
+  // ToDo This bail-out is non-deterministic, but the pseudo-clock
+  // model for INVERT can be used to make it deterministic
+  const bool check_for_timeout = this->time_limit_ < kHighsInf;
+  // To know when to bail out of buildKernel a model of how long
+  // buildKernel will take is needed
+  //
+  // If a pivot is deferred, nwork is increased to force another loop,
+  // so have to count the number of deferred pivots to know how many
+  // iterations might be needed
+  HighsInt num_defer_pivot = 0;
   // Initial timer frequency: may be reduced if iterations get slow
+  const HighsInt max_timer_frequency = 1000;
   HighsInt timer_frequency = 100;
-  double previous_iteration_time = 0;
+  // Need to maintain an averate iteration time
+  double previous_iteration_time = build_timer_->read();
   double average_iteration_time = 0;
-  const bool check_for_timeout = this->time_limit_ < kHighsInf;
-  HighsInt search_k = 0;
-
+  // Need to maintain an average rate of change of the number of
+  // active nonzeros - in order to predict when it will go to zero
+  HighsInt previous_num_active_nz = num_active_nz_;
+  double average_num_active_nz_change_rate = 0;
+  // Very occasionally, vectors of indices and values have to be moved
+  // in order to be resized, in which caes the cost of an iteration is
+  // abnormally high and needs to be excluded from the average
+  // calculation, so record when it happens
+  bool resize_data_shift = false;
+  // Parameters for the running average claculation
+  double mu0, mu1;
+
+  // Lambda functions that resize HighsInt/double vectors, with a
+  // check for them being moved as a result
+  auto resizeHighsInt = [&](std::vector<HighsInt>& i_vector,
+                            const HighsInt to_size) {
+    HighsInt* from_p = i_vector.data();
+    i_vector.resize(to_size);
+    if (i_vector.data() != from_p) resize_data_shift = true;
+  };
+
+  auto resizeDouble = [&](std::vector<double>& d_vector,
+                          const HighsInt to_size) {
+    double* from_p = d_vector.data();
+    d_vector.resize(to_size);
+    if (d_vector.data() != from_p) resize_data_shift = true;
+  };
+
+  // Work out the parameters for the running average claculation
+  auto runningAverageMu = [&]() {
+    mu0 = (1.0 * timer_frequency) / (10.0 * max_timer_frequency);
+    assert(mu0 <= 0.2);
+    mu1 = 1.0 - mu0;
+  };
+
+  runningAverageMu();
   const HighsInt check_nwork = -11;
   while (nwork-- > 0) {
     //    printf("\nnwork = %d\n", (int)nwork);
     if (nwork == check_nwork) {
       reportAsm();
     }
-    // Determine whether to return due to exceeding the time limit
-    if (check_for_timeout && search_k % timer_frequency == 0) {
+    // Determine whether to return due to (expecting to) exceed the
+    // time limit
+    if (check_for_timeout && search_k > 0 && search_k % timer_frequency == 0) {
+      // Get the current iteration time, and update the running
+      // average (if there has been no data shift on resize)
       double current_time = build_timer_->read();
       double time_difference = current_time - previous_iteration_time;
       previous_iteration_time = current_time;
       double iteration_time = time_difference / (1.0 * timer_frequency);
-      average_iteration_time =
-          0.9 * average_iteration_time + 0.1 * iteration_time;
-
-      if (time_difference > this->time_limit_ / 1e3)
-        timer_frequency = std::max(HighsInt(1), timer_frequency / 10);
-      HighsInt iterations_left = kernel_dim - search_k + 1;
+      if (!resize_data_shift) {
+        average_iteration_time =
+            mu1 * average_iteration_time + mu0 * iteration_time;
+        // Determine whether to reduce the timer frequency
+        if (timer_frequency > 1 && time_difference > this->time_limit_ / 1e1) {
+          timer_frequency = std::max(HighsInt(1), timer_frequency / 10);
+          runningAverageMu();
+        }
+      }
+      // Determine the current rate of change of the number of active
+      // nonzeros
+      double num_active_nz_change_rate =
+          static_cast<double>(num_active_nz_ - previous_num_active_nz) /
+          static_cast<double>(timer_frequency);
+      previous_num_active_nz = num_active_nz_;
+      average_num_active_nz_change_rate =
+          mu1 * average_num_active_nz_change_rate +
+          mu0 * num_active_nz_change_rate;
+
+      // Get an estimate of the number of iterations left, based on
+      // the bound and the average rate of change of the number of
+      // active nonzeros (if negative)
+      HighsInt iterations_left = kernel_dim - search_k + num_defer_pivot + 1;
+      if (average_num_active_nz_change_rate < -1) {
+        HighsInt active_nz_iterations_left =
+            -num_active_nz_ / average_num_active_nz_change_rate;
+        iterations_left = std::min(active_nz_iterations_left, iterations_left);
+      }
       double remaining_time_bound = average_iteration_time * iterations_left;
       double total_time_bound = current_time + remaining_time_bound;
+      // Update the record of bounds on total time - for logging in
+      // presolve
+      min_time_bound_ = std::min(total_time_bound, min_time_bound_);
+      max_time_bound_ = std::max(total_time_bound, max_time_bound_);
+      // Bail out if current or expected time exceeds the limit
       if (current_time > this->time_limit_ ||
           total_time_bound > this->time_limit_)
         return kBuildKernelReturnTimeout;
+      // Clear the record of vectors of indices and values having to
+      // be moved in order to be resized
+      resize_data_shift = false;
     }
-
     /**
      * 1. Search for the pivot
      */
@@ -1049,7 +1146,6 @@ HighsInt HFactor::buildKernel() {
                   (int)rank_deficiency);
       return rank_deficiency;
     }
-
     /**
      * 2. Elimination other elements by the pivot
      */
@@ -1061,6 +1157,9 @@ HighsInt HFactor::buildKernel() {
     //
     // Remove the pivot row index from the pivotal column of the
     // col-wise matrix. Also decreases the column count
+    //
+    // One active nonzero is lost
+    num_active_nz_--;
     double pivot_multiplier = colDelete(jColPivot, iRowPivot);
     // Remove the pivot column index from the pivotal row of the
     // row-wise matrix. Also decreases the row count
@@ -1078,6 +1177,7 @@ HighsInt HFactor::buildKernel() {
                   "Defer singular pivot = %11.4g\n", pivot_multiplier);
       // Matrix is singular, but defer return since other valid pivots
       // may exist.
+      num_defer_pivot++;
       assert(mr_count[iRowPivot] == original_pivotal_row_count - 1);
       if (mr_count[iRowPivot] == 0) {
         // The pivot corresponds to a singleton row. Entry is zeroed,
@@ -1123,6 +1223,8 @@ HighsInt HFactor::buildKernel() {
       mr_count_before[iRow] = mr_count[iRow];
       rowDelete(jColPivot, (int)iRow);
     }
+    // One active entry is lost for each entry in the pivotal column
+    num_active_nz_ -= (end_A - start_A);
     l_start.push_back(l_index.size());
     fake_fill += 2 * mc_count_a[jColPivot];
 
@@ -1149,6 +1251,8 @@ HighsInt HFactor::buildKernel() {
       const HighsInt my_end = my_start + my_count - 1;
       double my_pivot = colDelete(iCol, iRowPivot);
       colStoreN(iCol, iRowPivot, my_pivot);
+      // One active entry is lost for each entry in the pivotal row
+      num_active_nz_--;
 
       // 2.4.2. Elimination on the overlapping part
       HighsInt nFillin = mwz_column_count;
@@ -1167,6 +1271,8 @@ HighsInt HFactor::buildKernel() {
           mc_value[my_k] = value;
         }
       }
+      // One active entry is lost for each instance of cancellation
+      num_active_nz_ -= nCancel;
       fake_eliminate += mwz_column_count;
       fake_eliminate += nFillin * 2;
 
@@ -1196,8 +1302,8 @@ HighsInt HFactor::buildKernel() {
           mc_space[iCol] += max(mc_space[iCol], nFillin);
           HighsInt p5 = mc_start[iCol] = mc_index.size();
           HighsInt p7 = p5 + mc_space[iCol] - mc_count_n[iCol];
-          mc_index.resize(p5 + mc_space[iCol]);
-          mc_value.resize(p5 + mc_space[iCol]);
+          resizeHighsInt(mc_index, p5 + mc_space[iCol]);
+          resizeDouble(mc_value, p5 + mc_space[iCol]);
           copy(&mc_index[p1], &mc_index[p2], &mc_index[p5]);
           copy(&mc_value[p1], &mc_value[p2], &mc_value[p5]);
           copy(&mc_index[p3], &mc_index[p4], &mc_index[p7]);
@@ -1207,8 +1313,11 @@ HighsInt HFactor::buildKernel() {
         // 2.4.4.2 Fill into column copy
         for (HighsInt i = 0; i < mwz_column_count; i++) {
           HighsInt iRow = mwz_column_index[i];
-          if (mwz_column_mark[iRow])
+          if (mwz_column_mark[iRow]) {
             colInsert(iCol, iRow, -my_pivot * mwz_column_array[iRow]);
+            // One active entry is gained for each instance of fill-in
+            num_active_nz_++;
+          }
         }
 
         // 2.4.4.3 Fill into the row copy
@@ -1221,7 +1330,7 @@ HighsInt HFactor::buildKernel() {
               HighsInt p2 = p1 + mr_count[iRow];
               HighsInt p3 = mr_start[iRow] = mr_index.size();
               mr_space[iRow] *= 2;
-              mr_index.resize(p3 + mr_space[iRow]);
+              resizeHighsInt(mr_index, p3 + mr_space[iRow]);
               copy(&mr_index[p1], &mr_index[p2], &mr_index[p3]);
             }
             rowInsert(iCol, iRow);
@@ -1253,6 +1362,8 @@ HighsInt HFactor::buildKernel() {
         rlinkAdd(iRow, mr_count[iRow]);
       }
     }
+    // End of loop while(nwork-- > 0)
+    // Final execution has nwork = 0
   }
   build_synthetic_tick +=
       fake_search * 20 + fake_fill * 160 + fake_eliminate * 80;
diff --git a/highs/util/HFactor.h b/highs/util/HFactor.h
index 8e452846f8..3dfae4cfab 100644
--- a/highs/util/HFactor.h
+++ b/highs/util/HFactor.h
@@ -336,6 +336,9 @@ class HFactor {
   HighsInt invert_num_el;
   HighsInt kernel_dim;
   HighsInt kernel_num_el;
+  HighsInt num_active_nz_;
+  double min_time_bound_;
+  double max_time_bound_;
 
   /**
    * Data of the factor