From 59ba43306affb04e9e6dece7decbdc5bb54adfda Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Wed, 25 Mar 2026 12:17:43 +0000 Subject: [PATCH 01/16] Added clocks to presolve analysis --- cmake/sources-python.cmake | 1 + cmake/sources.cmake | 1 + highs/lp_data/HConst.h | 4 +++- highs/mip/MipTimer.h | 2 +- highs/presolve/HPresolve.cpp | 2 +- highs/presolve/HPresolveAnalysis.cpp | 14 +++++++++++++- highs/presolve/HPresolveAnalysis.h | 20 ++++++++++++++++---- 7 files changed, 36 insertions(+), 8 deletions(-) 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 555166a026..ddbdedbb6a 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -546,6 +546,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/lp_data/HConst.h b/highs/lp_data/HConst.h index 020a91504a..67071baa59 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..39c89bc14f 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, diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index c27404da9e..3bc3c1a222 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -5746,7 +5746,7 @@ 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); if (options->presolve != kHighsOffString) { if (mipsolver) mipsolver->mipdata_->cliquetable.setPresolveFlag(true); diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index f4724ebd05..c452e21ea3 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -7,17 +7,29 @@ /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ #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; + 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) { diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h index 9fa48f32d2..6960abee3e 100644 --- a/highs/presolve/HPresolveAnalysis.h +++ b/highs/presolve/HPresolveAnalysis.h @@ -8,10 +8,18 @@ /**@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; @@ -32,6 +40,9 @@ class HPresolveAnalysis { HighsInt num_deleted_rows0_; HighsInt num_deleted_cols0_; HighsPresolveLog presolve_log_; + + HighsTimerClock presolve_clocks_; + bool analyse_presolve_time_; // for LP presolve // @@ -39,7 +50,8 @@ class HPresolveAnalysis { // 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); @@ -49,4 +61,4 @@ class HPresolveAnalysis { friend class HPresolve; }; -#endif +#endif /* PRESOLVE_HPRESOLVEANALYSIS_H_ */ From 697a00ee6e9bb98616a33eaf0a45bf20515023c5 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Wed, 25 Mar 2026 12:43:23 +0000 Subject: [PATCH 02/16] Now to see why top presolve clock isn't recording/reporting --- highs/mip/MipTimer.h | 26 +++++++++++++------------- highs/presolve/HPresolve.cpp | 5 ++++- highs/presolve/HPresolveAnalysis.cpp | 19 +++++++++++++++++++ highs/presolve/HPresolveAnalysis.h | 3 +++ 4 files changed, 39 insertions(+), 14 deletions(-) diff --git a/highs/mip/MipTimer.h b/highs/mip/MipTimer.h index 39c89bc14f..1cb06d0068 100644 --- a/highs/mip/MipTimer.h +++ b/highs/mip/MipTimer.h @@ -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 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 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) { @@ -443,7 +443,7 @@ class MipTimer { kMipClockRootSeparationCentralRounding, kMipClockRootSeparationEvaluateRootLp}; reportMipClockList("MipRootSeparation", mip_clock_list, mip_timer_clock, - kMipClockRootSeparation); //, tolerance_percent_report); + kMipClockRootSeparation); //, kMipClockTolerancePercentReport); }; void reportMipSearchClock(const HighsTimerClock& mip_timer_clock) { @@ -455,7 +455,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 +463,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 +472,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) { @@ -482,7 +482,7 @@ class MipTimer { kMipClockOpenNodesToQueue1, kMipClockEvaluateNode1, kMipClockNodeSearchSeparation}; //, kMipClockStoreBasis}; reportMipClockList("MipNodeSearch", mip_clock_list, mip_timer_clock, - kMipClockNodeSearch); //, tolerance_percent_report); + kMipClockNodeSearch); //, kMipClockTolerancePercentReport); }; void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) { @@ -490,7 +490,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/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 3bc3c1a222..8edfa6e518 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" @@ -5747,7 +5748,7 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { // effectiveness analysis_.setup(this->model, this->options, this->numDeletedRows, this->numDeletedCols, silent, this->timer); - + analysis_.presolveTimerStart(kPresolveClockPresolve); if (options->presolve != kHighsOffString) { if (mipsolver) mipsolver->mipdata_->cliquetable.setPresolveFlag(true); @@ -5953,6 +5954,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 diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index c452e21ea3..a96a4b0d5e 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -249,3 +249,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(presolve_clocks_); +} + diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h index 6960abee3e..b6ba39fb00 100644 --- a/highs/presolve/HPresolveAnalysis.h +++ b/highs/presolve/HPresolveAnalysis.h @@ -58,6 +58,9 @@ class HPresolveAnalysis { 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; }; From 2d8f4b2b6587bde999959b25bcaf1cf03d696590 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Wed, 25 Mar 2026 12:50:11 +0000 Subject: [PATCH 03/16] Added highs/presolve/PresolveTimer.h and unit test presolve-time --- check/TestPresolve.cpp | 10 ++++++ highs/presolve/PresolveTimer.h | 64 ++++++++++++++++++++++++++++++++++ 2 files changed, 74 insertions(+) create mode 100644 highs/presolve/PresolveTimer.h diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index 8af6618f34..9820494a04 100644 --- a/check/TestPresolve.cpp +++ b/check/TestPresolve.cpp @@ -910,3 +910,13 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") { REQUIRE(highs.presolve() == HighsStatus::kOk); REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible); } + +TEST_CASE("presolve-time", "[highs_test_presolve]") { + std::string model_file = + std::string(HIGHS_DIR) + "/check/instances/shell.mps"; + Highs highs; + // highs.setOptionValue("output_flag", dev_run); + highs.readModel(model_file); + highs.setOptionValue("highs_analysis_level", 256); + REQUIRE(highs.presolve() == HighsStatus::kOk); +} diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h new file mode 100644 index 0000000000..23622b54a4 --- /dev/null +++ b/highs/presolve/PresolveTimer.h @@ -0,0 +1,64 @@ +/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ +/* */ +/* 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, + 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& clock = presolve_timer_clock.clock_; + + clock.resize(kNumPresolveClock); + clock[kPresolveClockTotal] = 0; + clock[kPresolveClockPresolve] = timer_pointer->clock_def("Presolve"); + }; + + bool reportPresolveClockList(const char* grepStamp, + const std::vector presolve_clock_list, + const HighsTimerClock& presolve_timer_clock, + const HighsInt kPresolveClockIdeal = kPresolveClockTotal, + const double tolerance_percent_report_ = -1) { + HighsTimer* timer_pointer = presolve_timer_clock.timer_pointer_; + if (!timer_pointer->printf_flag) return false; + const std::vector& clock = presolve_timer_clock.clock_; + HighsInt presolve_clock_list_size = presolve_clock_list.size(); + std::vector 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 reportPresolveCoreClock(const HighsTimerClock& presolve_timer_clock) { + const std::vector presolve_clock_list{ + kPresolveClockPresolve}; + reportPresolveClockList("PresolveCore_", presolve_clock_list, presolve_timer_clock, + kPresolveClockTotal); + }; +}; + +#endif /* PRESOLVE_PRESOLVETIMER_H_ */ From d9862d4da4ba0f597b9952b3ad76d836361233a8 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Wed, 25 Mar 2026 22:20:27 +0000 Subject: [PATCH 04/16] Logging temporary data in dependent equations reduction --- check/TestPdlpHi.cpp | 1 - check/TestPresolve.cpp | 2 +- highs/mip/MipTimer.h | 10 +- highs/pdlp/CupdlpWrapper.cpp | 7 +- highs/pdlp/hipdlp/defs.hpp | 14 +- highs/pdlp/hipdlp/logger.cc | 21 +- highs/pdlp/hipdlp/logger.hpp | 14 +- highs/pdlp/hipdlp/pdhg.cc | 540 +++++++++++++++------------ highs/pdlp/hipdlp/pdhg.hpp | 132 ++++--- highs/pdlp/hipdlp/scaling.cc | 12 +- highs/pdlp/hipdlp/scaling.hpp | 5 +- highs/presolve/HPresolve.cpp | 86 ++++- highs/presolve/HPresolveAnalysis.cpp | 12 +- highs/presolve/HPresolveAnalysis.h | 8 +- highs/presolve/PresolveTimer.h | 117 +++++- 15 files changed, 622 insertions(+), 359 deletions(-) diff --git a/check/TestPdlpHi.cpp b/check/TestPdlpHi.cpp index 1841f1a548..27d415d3de 100644 --- a/check/TestPdlpHi.cpp +++ b/check/TestPdlpHi.cpp @@ -38,4 +38,3 @@ TEST_CASE("hi-pdlp", "[pdlp]") { } h.resetGlobalScheduler(true); } - diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index 9820494a04..f8478a0fe2 100644 --- a/check/TestPresolve.cpp +++ b/check/TestPresolve.cpp @@ -912,7 +912,7 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") { } TEST_CASE("presolve-time", "[highs_test_presolve]") { - std::string model_file = + std::string model_file = std::string(HIGHS_DIR) + "/check/instances/shell.mps"; Highs highs; // highs.setOptionValue("output_flag", dev_run); diff --git a/highs/mip/MipTimer.h b/highs/mip/MipTimer.h index 1cb06d0068..3d55feec73 100644 --- a/highs/mip/MipTimer.h +++ b/highs/mip/MipTimer.h @@ -442,8 +442,9 @@ class MipTimer { kMipClockRootSeparationFinishAnalyticCentreComputation, kMipClockRootSeparationCentralRounding, kMipClockRootSeparationEvaluateRootLp}; - reportMipClockList("MipRootSeparation", mip_clock_list, mip_timer_clock, - kMipClockRootSeparation); //, kMipClockTolerancePercentReport); + reportMipClockList( + "MipRootSeparation", mip_clock_list, mip_timer_clock, + kMipClockRootSeparation); //, kMipClockTolerancePercentReport); }; void reportMipSearchClock(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); //, kMipClockTolerancePercentReport); + reportMipClockList( + "MipNodeSearch", mip_clock_list, mip_timer_clock, + kMipClockNodeSearch); //, kMipClockTolerancePercentReport); }; void reportMipSeparationClock(const HighsTimerClock& mip_timer_clock) { diff --git a/highs/pdlp/CupdlpWrapper.cpp b/highs/pdlp/CupdlpWrapper.cpp index f59b0b2f3e..d9cb5158a6 100644 --- a/highs/pdlp/CupdlpWrapper.cpp +++ b/highs/pdlp/CupdlpWrapper.cpp @@ -211,13 +211,12 @@ HighsStatus solveLpCupdlp(const HighsOptions& options, HighsTimer& timer, highs_solution.col_value.data(), lp.num_col_, highs_solution.row_dual.data(), lp.num_row_); if (w->debug_pdlp_log_file_) fclose(w->debug_pdlp_log_file_); - // Moved this from LP_SolvePDHG so w->debug_pdlp_log_file_ can - // still be used - // + // Moved this from LP_SolvePDHG so w->debug_pdlp_log_file_ can + // still be used + // #endif PDHG_Destroy(&w); - model_status = HighsModelStatus::kUnknown; highs_solution.value_valid = value_valid; highs_solution.dual_valid = dual_valid; diff --git a/highs/pdlp/hipdlp/defs.hpp b/highs/pdlp/hipdlp/defs.hpp index 1d740d3658..3b97f598fa 100644 --- a/highs/pdlp/hipdlp/defs.hpp +++ b/highs/pdlp/hipdlp/defs.hpp @@ -20,7 +20,12 @@ enum class Device { CPU, GPU }; enum class ScalingMethod { NONE, RUIZ, POCK_CHAMBOLLE, L2_NORM, COMBINED }; -enum class RestartStrategy { NO_RESTART, FIXED_RESTART, ADAPTIVE_RESTART, HALPERN_RESTART}; +enum class RestartStrategy { + NO_RESTART, + FIXED_RESTART, + ADAPTIVE_RESTART, + HALPERN_RESTART +}; enum class StepSizeStrategy { FIXED, ADAPTIVE, MALITSKY_POCK, PID }; @@ -67,7 +72,8 @@ struct PrimalDualParams { HighsInt fixed_restart_interval; bool use_halpern_restart = false; - double halpern_gamma = 1.0; // 0: standard Halpern, 1: full reflection, in between: over relaxation + double halpern_gamma = 1.0; // 0: standard Halpern, 1: full reflection, in + // between: over relaxation // Scaling parameters bool use_ruiz_scaling = false; @@ -191,8 +197,8 @@ struct DetailedTimings { " - Matrix multiply: %6.2f s (%3.0f)\n", matrix_multiply_time, matrix_multiply_time / total_time * 100); highsLogUser(log_options, HighsLogType::kInfo, - " - Projection: %6.2f s (%3.0f)\n", - projection_time, projection_time / total_time * 100); + " - Projection: %6.2f s (%3.0f)\n", projection_time, + projection_time / total_time * 100); highsLogUser(log_options, HighsLogType::kInfo, " - Step size adjust: %6.2f s (%3.0f)\n", step_size_adjustment_time, diff --git a/highs/pdlp/hipdlp/logger.cc b/highs/pdlp/hipdlp/logger.cc index a54d3b0cb4..4578d3ccc5 100644 --- a/highs/pdlp/hipdlp/logger.cc +++ b/highs/pdlp/hipdlp/logger.cc @@ -40,8 +40,12 @@ void Logger::log(LogLevel level, const std::string& message) const { highsLogUser(log_options_, HighsLogType::kInfo, "%s\n", message.c_str()); } -void Logger::info(const std::string& message) const { log(LogLevel::kInfo, message); } -void Logger::detailed(const std::string& message) const { log(LogLevel::kDetailed, message); } +void Logger::info(const std::string& message) const { + log(LogLevel::kInfo, message); +} +void Logger::detailed(const std::string& message) const { + log(LogLevel::kDetailed, message); +} void Logger::verbose(const std::string& message) const { log(LogLevel::kVerbose, message); } @@ -60,8 +64,7 @@ void Logger::print_params(const PrimalDualParams& params) const { std::map restart_map = { {RestartStrategy::NO_RESTART, "None"}, {RestartStrategy::FIXED_RESTART, "Fixed"}, - {RestartStrategy::ADAPTIVE_RESTART, "Adaptive"} - }; + {RestartStrategy::ADAPTIVE_RESTART, "Adaptive"}}; std::map step_size_map = { {StepSizeStrategy::FIXED, "Fixed"}, {StepSizeStrategy::ADAPTIVE, "Adaptive"}, @@ -120,14 +123,16 @@ void Logger::print_params(const PrimalDualParams& params) const { } void Logger::print_iteration_header() const { - info(" Iter Primal Feas Dual Feas P-D Gap Step Size Time"); + info( + " Iter Primal Feas Dual Feas P-D Gap Step Size " + "Time"); } void Logger::print_iteration_stats(HighsInt iter, const SolverResults& results, const double step_size, - const double time) const { + const double time) const { std::stringstream ss; -// clang-format off + // clang-format off ss << std::fixed << std::setprecision(9) << std::setw(9) << iter << " " << std::scientific << std::setprecision(2) << std::setw(11) << results.primal_feasibility << " " @@ -137,7 +142,7 @@ void Logger::print_iteration_stats(HighsInt iter, const SolverResults& results, << std::setw(9) << step_size << " " << std::setprecision(1) << std::setw(8) << time; -// clang-format on + // clang-format on info(ss.str()); } diff --git a/highs/pdlp/hipdlp/logger.hpp b/highs/pdlp/hipdlp/logger.hpp index c943f9bee1..99372034a7 100644 --- a/highs/pdlp/hipdlp/logger.hpp +++ b/highs/pdlp/hipdlp/logger.hpp @@ -23,11 +23,11 @@ // Log verbosity level enum class LogLevel { - kNone, // No output - kInfo, // Standard output: summary, termination, major events - kDetailed, // Developer output - kVerbose, // Detailed output: iteration-level info - kDebug // Verbose + debug info for developers + kNone, // No output + kInfo, // Standard output: summary, termination, major events + kDetailed, // Developer output + kVerbose, // Detailed output: iteration-level info + kDebug // Verbose + debug info for developers }; class Logger { @@ -48,7 +48,8 @@ class Logger { void print_header() const; void print_params(const PrimalDualParams& params) const; void print_iteration_header() const; - void print_iteration_stats(HighsInt iter, const SolverResults& current_results, + void print_iteration_stats(HighsInt iter, + const SolverResults& current_results, const double current_eta, const double time) const; void print_summary(const SolverResults& results, HighsInt total_iter, double total_time) const; @@ -56,6 +57,7 @@ class Logger { log_options_ = log_options; } LogLevel getConsoleLevel() const { return console_level_; } + private: void log(LogLevel level, const std::string& message) const; LogLevel console_level_; diff --git a/highs/pdlp/hipdlp/pdhg.cc b/highs/pdlp/hipdlp/pdhg.cc index fd8981ed06..1f4cd604ae 100644 --- a/highs/pdlp/hipdlp/pdhg.cc +++ b/highs/pdlp/hipdlp/pdhg.cc @@ -53,7 +53,7 @@ void PDLPSolver::printConstraintInfo() { // Count constraint types BEFORE preprocessing HighsInt eq_count = 0, leq_count = 0, geq_count = 0, bound_count = 0, - free_count = 0; + free_count = 0; for (HighsInt i = 0; i < nRows; ++i) { bool has_lower = original_lp_->row_lower_[i] > -kHighsInf; @@ -99,7 +99,7 @@ void PDLPSolver::printConstraintInfo() { logger_.detailed("=== BEFORE PREPROCESSING ==="); logger_.detailed("Rows: " + std::to_string(nRows) + - ", Cols: " + std::to_string(nCols)); + ", Cols: " + std::to_string(nCols)); logger_.detailed("\nConstraint types:"); logger_.detailed(" Equality constraints (=): " + std::to_string(eq_count)); logger_.detailed(" One-sided inequality (>=): " + std::to_string(geq_count)); @@ -109,22 +109,23 @@ void PDLPSolver::printConstraintInfo() { logger_.detailed("\nVariable bound types:"); logger_.detailed(" Fixed variables (l = u): " + std::to_string(var_fixed)); - logger_.detailed(" Boxed variables (l <= x <= u): " + std::to_string(var_boxed)); + logger_.detailed(" Boxed variables (l <= x <= u): " + + std::to_string(var_boxed)); logger_.detailed(" Lower bounded only (l <= x): " + - std::to_string(var_lower_only)); + std::to_string(var_lower_only)); logger_.detailed(" Upper bounded only (x <= u): " + - std::to_string(var_upper_only)); + std::to_string(var_upper_only)); logger_.detailed(" Free variables: " + std::to_string(var_free)); logger_.detailed("\n=== AFTER PREPROCESSING ==="); logger_.detailed("Rows: " + std::to_string(lp_.num_row_) + - ", Cols: " + std::to_string(lp_.num_col_)); + ", Cols: " + std::to_string(lp_.num_col_)); logger_.detailed("Equality rows (first " + std::to_string(num_eq_rows_) + - " rows)"); + " rows)"); logger_.detailed("Inequality rows (remaining " + - std::to_string(lp_.num_row_ - num_eq_rows_) + " rows)"); + std::to_string(lp_.num_row_ - num_eq_rows_) + " rows)"); logger_.detailed("Slack variables added: " + - std::to_string(lp_.num_col_ - nCols)); + std::to_string(lp_.num_col_ - nCols)); // Count variable bounds in processed LP HighsInt proc_var_fixed = 0, proc_var_lower_only = 0, proc_var_upper_only = 0; @@ -152,8 +153,10 @@ void PDLPSolver::printConstraintInfo() { logger_.detailed("\nProcessed variable bound types:"); logger_.detailed(" Fixed variables: " + std::to_string(proc_var_fixed)); logger_.detailed(" Boxed variables: " + std::to_string(proc_var_boxed)); - logger_.detailed(" Lower bounded only: " + std::to_string(proc_var_lower_only)); - logger_.detailed(" Upper bounded only: " + std::to_string(proc_var_upper_only)); + logger_.detailed(" Lower bounded only: " + + std::to_string(proc_var_lower_only)); + logger_.detailed(" Upper bounded only: " + + std::to_string(proc_var_upper_only)); logger_.detailed(" Free variables: " + std::to_string(proc_var_free)); } @@ -352,10 +355,11 @@ void PDLPSolver::preprocessLp() { unscaled_rhs_norm_ = linalg::vector_norm(processed_lp.row_lower_); logger_.detailed("Preprocessing complete. New dimensions: " + - std::to_string(processed_lp.num_row_) + " rows, " + - std::to_string(processed_lp.num_col_) + " cols."); - logger_.detailed("Unscaled norms: ||c|| = " + std::to_string(unscaled_c_norm_) + - ", ||b|| = " + std::to_string(unscaled_rhs_norm_)); + std::to_string(processed_lp.num_row_) + " rows, " + + std::to_string(processed_lp.num_col_) + " cols."); + logger_.detailed( + "Unscaled norms: ||c|| = " + std::to_string(unscaled_c_norm_) + + ", ||b|| = " + std::to_string(unscaled_rhs_norm_)); printConstraintInfo(); #if PDLP_PROFILE @@ -461,8 +465,8 @@ PostSolveRetcode PDLPSolver::postprocess(HighsSolution& solution) { for (HighsInt col = 0; col < original_num_col_; ++col) { double x_val = x_current_[col]; // Use unscaled x values - for (HighsInt el = orig_matrix.start_[col]; el < orig_matrix.start_[col + 1]; - ++el) { + for (HighsInt el = orig_matrix.start_[col]; + el < orig_matrix.start_[col + 1]; ++el) { HighsInt row = orig_matrix.index_[el]; double a_val = orig_matrix.value_[el]; ax_original[row] += a_val * x_val; @@ -502,20 +506,21 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { #if PDLP_PROFILE hipdlpTimerStart(kHipdlpClockSolve); #endif - + // 1. Initialization and Setup #ifdef CUPDLP_GPU - setupGpu(); + setupGpu(); #endif initializeStepSizes(); - initialize(); // Resets vectors, caches, and sets initial x_current, y_current + initialize(); // Resets vectors, caches, and sets initial x_current, + // y_current best_primal_weight_ = primal_weight_; best_primal_dual_residual_gap_ = std::numeric_limits::infinity(); // Initialize internal solver state restart_scheme_.passParams(¶ms_); restart_scheme_.Initialize(results_); - + // Copy initial input x,y to internal state x_current_ = x; y_current_ = y; @@ -525,8 +530,12 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { y_anchor_ = y_current_; halpern_iteration_ = 0; #ifdef CUPDLP_GPU - CUDA_CHECK(cudaMemcpy(d_x_anchor_, d_x_current_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_y_anchor_, d_y_current_, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_x_anchor_, d_x_current_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_y_anchor_, d_y_current_, + lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToDevice)); #endif } @@ -537,8 +546,10 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { CUDA_CHECK(cudaMemset(d_x_avg_, 0, lp_.num_col_ * sizeof(double))); CUDA_CHECK(cudaMemset(d_y_avg_, 0, lp_.num_row_ * sizeof(double))); sum_weights_gpu_ = 0.0; - CUDA_CHECK(cudaMemcpy(d_x_current_, x.data(), lp_.num_col_ * sizeof(double), cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_y_current_, y.data(), lp_.num_row_ * sizeof(double), cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_x_current_, x.data(), lp_.num_col_ * sizeof(double), + cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_y_current_, y.data(), lp_.num_row_ * sizeof(double), + cudaMemcpyHostToDevice)); linalgGpuAx(d_x_current_, d_ax_current_); linalgGpuATy(d_y_current_, d_aty_current_); #else @@ -548,14 +559,14 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { #endif logger_.print_iteration_header(); - + // termination_status is not known. Setting it to NOTSET rather than // OPTIMAL means that if it's not set elsewhere then optimality is // not returned in error TerminationStatus termination_status = TerminationStatus::NOTSET; bool do_restart = false; double last_trial_fpe = std::numeric_limits::infinity(); - + final_iter_count_ = 0; #ifdef CUPDLP_GPU @@ -574,13 +585,15 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { // -- Step 1 (Major, isolated for restart-FPE check) -- #ifdef CUPDLP_GPU - CUDA_CHECK(cudaMemcpyAsync(d_primal_step_size_, &stepsize_.primal_step, - sizeof(double), cudaMemcpyHostToDevice, - gpu_stream_)); - CUDA_CHECK(cudaMemcpyAsync(d_dual_step_size_, &stepsize_.dual_step, - sizeof(double), cudaMemcpyHostToDevice, - gpu_stream_)); - CUDA_CHECK(cudaMemcpyAsync(d_halpern_iteration_, &halpern_iteration_, sizeof(int), cudaMemcpyHostToDevice, gpu_stream_)); + CUDA_CHECK(cudaMemcpyAsync(d_primal_step_size_, &stepsize_.primal_step, + sizeof(double), cudaMemcpyHostToDevice, + gpu_stream_)); + CUDA_CHECK(cudaMemcpyAsync(d_dual_step_size_, &stepsize_.dual_step, + sizeof(double), cudaMemcpyHostToDevice, + gpu_stream_)); + CUDA_CHECK(cudaMemcpyAsync(d_halpern_iteration_, &halpern_iteration_, + sizeof(int), cudaMemcpyHostToDevice, + gpu_stream_)); performHalpernPdhgStepGpu(true, 1); #else performHalpernPdhgStep(true, 1); @@ -593,17 +606,23 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { fpe_ = computeFixedPointError(); #endif initial_fpe_ = fpe_; - // if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << final_iter_count_ << " | Restarted with initial_fpe_ = " << std::setprecision(10) << initial_fpe_ << std::endl; + // if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << final_iter_count_ << " + // | Restarted with initial_fpe_ = " << std::setprecision(10) << + // initial_fpe_ << std::endl; do_restart = false; } // copy back x from gpu and print out #ifdef CUPDLP_GPU - if(final_iter_count_ == 9) { + if (final_iter_count_ == 9) { std::vector x_next_host(lp_.num_col_); - CUDA_CHECK(cudaMemcpy(x_next_host.data(), d_pdhg_primal_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(x_next_host.data(), d_pdhg_primal_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); std::cout << "x_next (before 9): "; - for(int i=0; i<3; i++) std::cout << std::fixed << std::setprecision(10) << x_next_host[i] << " "; + for (int i = 0; i < 3; i++) + std::cout << std::fixed << std::setprecision(10) << x_next_host[i] + << " "; std::cout << "\n"; } #endif @@ -635,41 +654,52 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { #endif #ifdef CUPDLP_GPU - if(final_iter_count_ == 9) { + if (final_iter_count_ == 9) { std::vector x_next_host(lp_.num_col_); std::vector x_current_host(lp_.num_col_); - CUDA_CHECK(cudaMemcpy(x_current_host.data(), d_x_current_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); - CUDA_CHECK(cudaMemcpy(x_next_host.data(), d_pdhg_primal_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(x_current_host.data(), d_x_current_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(x_next_host.data(), d_pdhg_primal_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); std::cout << "x_next (after 9): "; - for(int i=0; i<3; i++) std::cout << std::fixed << std::setprecision(10) << x_next_host[i] << " "; + for (int i = 0; i < 3; i++) + std::cout << std::fixed << std::setprecision(10) << x_next_host[i] + << " "; std::cout << "\n"; std::cout << "x_current (after 9): "; - for(int i=0; i<3; i++) std::cout << std::fixed << std::setprecision(10) << x_current_host[i] << " "; + for (int i = 0; i < 3; i++) + std::cout << std::fixed << std::setprecision(10) << x_current_host[i] + << " "; std::cout << "\n"; } #endif - // Compute Error for Restart Check - #ifdef CUPDLP_GPU - fpe_ = computeFixedPointErrorGpu(); - #else - fpe_ = computeFixedPointError(); - #endif -//if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << final_iter_count_ << " | FPE at check: " << std::setprecision(10) << fpe_ << std::endl; - +// Compute Error for Restart Check +#ifdef CUPDLP_GPU + fpe_ = computeFixedPointErrorGpu(); +#else + fpe_ = computeFixedPointError(); +#endif + // if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << final_iter_count_ << " + // | FPE at check: " << std::setprecision(10) << fpe_ << std::endl; + halpern_iteration_ += PDHG_CHECK_INTERVAL; final_iter_count_ += PDHG_CHECK_INTERVAL; // -- Convergence Check -- TerminationStatus check_status; - bool should_terminate = runConvergenceCheckAndRestart(final_iter_count_, x, y, check_status); - + bool should_terminate = + runConvergenceCheckAndRestart(final_iter_count_, x, y, check_status); + if (should_terminate) { termination_status = check_status; - break; + break; } - // -- Check Adaptive Restart (Mimicking `should_do_adaptive_restart` in solver.cu) -- + // -- Check Adaptive Restart (Mimicking `should_do_adaptive_restart` in + // solver.cu) -- if (final_iter_count_ == PDHG_CHECK_INTERVAL) { do_restart = true; } else if (final_iter_count_ > PDHG_CHECK_INTERVAL) { @@ -682,8 +712,7 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { } } if (halpern_iteration_ >= - restart_scheme_.GetArtificialRestartThreshold() * - final_iter_count_) { + restart_scheme_.GetArtificialRestartThreshold() * final_iter_count_) { do_restart = true; } } @@ -693,13 +722,23 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { if (do_restart) { if (params_.step_size_strategy == StepSizeStrategy::PID) { updatePrimalWeightAtRestart(results_); -//if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << final_iter_count_ << " | Updated primal_weight_ to " << primal_weight_ << std::endl; + // if (DEBUG_MODE) std::cout << "[DEBUG] iter=" << + // final_iter_count_ << " | Updated primal_weight_ to " << + // primal_weight_ << std::endl; } #ifdef CUPDLP_GPU - CUDA_CHECK(cudaMemcpy(d_x_anchor_, d_pdhg_primal_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_y_anchor_, d_pdhg_dual_, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_x_current_, d_pdhg_primal_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_y_current_, d_pdhg_dual_, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_x_anchor_, d_pdhg_primal_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_y_anchor_, d_pdhg_dual_, + lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_x_current_, d_pdhg_primal_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_y_current_, d_pdhg_dual_, + lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToDevice)); #else x_anchor_ = x_next_; y_anchor_ = y_next_; @@ -708,7 +747,7 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { linalg::Ax(lp_, x_current_, Ax_cache_); linalg::ATy(lp_, y_current_, ATy_cache_); #endif - + halpern_iteration_ = 0; last_trial_fpe = std::numeric_limits::infinity(); } @@ -725,11 +764,11 @@ void PDLPSolver::solve(std::vector& x, std::vector& y) { logger_.info("Iteration limit reached without convergence"); termination_status = TerminationStatus::MAXITER; } - + solveReturn(termination_status); } -double PDLPSolver::computeFixedPointError(){ +double PDLPSolver::computeFixedPointError() { double primal_norm_sq = 0.0; double dual_norm_sq = 0.0; double cross_term = 0.0; @@ -754,7 +793,8 @@ double PDLPSolver::computeFixedPointError(){ cross_term += delta_x[i] * AT_delta_y[i]; } - double movement = primal_norm_sq * params_.omega + dual_norm_sq / params_.omega; + double movement = + primal_norm_sq * params_.omega + dual_norm_sq / params_.omega; double interaction = 2.0 * params_.eta * cross_term; return std::sqrt(std::max(0.0, movement + interaction)); @@ -763,45 +803,55 @@ double PDLPSolver::computeFixedPointError(){ #ifdef CUPDLP_GPU double PDLPSolver::computeFixedPointErrorGpu() { double alpha_minus_one = -1.0; - - // 1. delta_x = x_next_ - reflected_x_ - // (Assuming d_pdhg_primal_ maps to x_next_ and d_x_next_ is used as reflected_x_ in your minor/major steps) - CUDA_CHECK(cudaMemcpy(d_delta_x_, d_pdhg_primal_, a_num_cols_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUBLAS_CHECK(cublasDaxpy(cublas_handle_, a_num_cols_, &alpha_minus_one, d_x_next_, 1, d_delta_x_, 1)); + + // 1. delta_x = x_next_ - reflected_x_ + // (Assuming d_pdhg_primal_ maps to x_next_ and d_x_next_ is used as + // reflected_x_ in your minor/major steps) + CUDA_CHECK(cudaMemcpy(d_delta_x_, d_pdhg_primal_, + a_num_cols_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUBLAS_CHECK(cublasDaxpy(cublas_handle_, a_num_cols_, &alpha_minus_one, + d_x_next_, 1, d_delta_x_, 1)); // 2. delta_y = y_next_ - reflected_y_ - CUDA_CHECK(cudaMemcpy(d_delta_y_, d_pdhg_dual_, a_num_rows_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUBLAS_CHECK(cublasDaxpy(cublas_handle_, a_num_rows_, &alpha_minus_one, d_y_next_, 1, d_delta_y_, 1)); + CUDA_CHECK(cudaMemcpy(d_delta_y_, d_pdhg_dual_, a_num_rows_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUBLAS_CHECK(cublasDaxpy(cublas_handle_, a_num_rows_, &alpha_minus_one, + d_y_next_, 1, d_delta_y_, 1)); // 3. AT_delta_y = A^T * delta_y linalgGpuATy(d_delta_y_, d_AT_delta_y_); // 4. Compute norms and cross term double primal_norm = 0.0, dual_norm = 0.0, cross_term = 0.0; - - CUBLAS_CHECK(cublasDnrm2(cublas_handle_, a_num_cols_, d_delta_x_, 1, &primal_norm)); - CUBLAS_CHECK(cublasDnrm2(cublas_handle_, a_num_rows_, d_delta_y_, 1, &dual_norm)); - CUBLAS_CHECK(cublasDdot(cublas_handle_, a_num_cols_, d_delta_x_, 1, d_AT_delta_y_, 1, &cross_term)); + + CUBLAS_CHECK( + cublasDnrm2(cublas_handle_, a_num_cols_, d_delta_x_, 1, &primal_norm)); + CUBLAS_CHECK( + cublasDnrm2(cublas_handle_, a_num_rows_, d_delta_y_, 1, &dual_norm)); + CUBLAS_CHECK(cublasDdot(cublas_handle_, a_num_cols_, d_delta_x_, 1, + d_AT_delta_y_, 1, &cross_term)); double primal_norm_sq = primal_norm * primal_norm; double dual_norm_sq = dual_norm * dual_norm; - double movement = primal_norm_sq * params_.omega + dual_norm_sq / params_.omega; + double movement = + primal_norm_sq * params_.omega + dual_norm_sq / params_.omega; double interaction = 2.0 * params_.eta * cross_term; return std::sqrt(std::max(0.0, movement + interaction)); } #endif -bool PDLPSolver::runConvergenceCheckAndRestart(size_t iter, - std::vector& output_x, +bool PDLPSolver::runConvergenceCheckAndRestart(size_t iter, + std::vector& output_x, std::vector& output_y, TerminationStatus& status) { // 1. Compute Average Iterate (GPU or CPU) #ifdef CUPDLP_GPU computeAverageIterateGpu(); #else - computeAverageIterate(Ax_avg_, ATy_avg_); + computeAverageIterate(Ax_avg_, ATy_avg_); #endif SolverResults current_results; @@ -812,28 +862,25 @@ bool PDLPSolver::runConvergenceCheckAndRestart(size_t iter, // 2. Check Convergence (GPU or CPU) #ifdef CUPDLP_GPU if (params_.use_halpern_restart && d_pdhg_primal_ != nullptr) { - // For Halpern, "current iterate" convergence uses the PDHG projected iterates - // Need Ax and ATy of pdhg iterates: - linalgGpuAx(d_pdhg_primal_, d_ax_next_); // reuse d_ax_next_ as scratch - linalgGpuATy(d_pdhg_dual_, d_aty_next_); // reuse d_aty_next_ as scratch + // For Halpern, "current iterate" convergence uses the PDHG projected + // iterates Need Ax and ATy of pdhg iterates: + linalgGpuAx(d_pdhg_primal_, d_ax_next_); // reuse d_ax_next_ as scratch + linalgGpuATy(d_pdhg_dual_, d_aty_next_); // reuse d_aty_next_ as scratch current_converged = checkConvergenceGpu( - iter, d_pdhg_primal_, d_pdhg_dual_, - d_ax_next_, d_aty_next_, - params_.tolerance, current_results, "[L-GPU]", - d_dSlackPos_, d_dSlackNeg_, - true); // <-- TRUE: Use the Halpern slack saved in d_dSlackPos_ + iter, d_pdhg_primal_, d_pdhg_dual_, d_ax_next_, d_aty_next_, + params_.tolerance, current_results, "[L-GPU]", d_dSlackPos_, + d_dSlackNeg_, + true); // <-- TRUE: Use the Halpern slack saved in d_dSlackPos_ } else { - current_converged = checkConvergenceGpu( - iter, d_x_current_, d_y_current_, - d_ax_current_, d_aty_current_, - params_.tolerance, current_results, "[L-GPU]", - d_dSlackPos_, d_dSlackNeg_, - false); // <-- FALSE: Recompute standard slack + current_converged = + checkConvergenceGpu(iter, d_x_current_, d_y_current_, d_ax_current_, + d_aty_current_, params_.tolerance, current_results, + "[L-GPU]", d_dSlackPos_, d_dSlackNeg_, + false); // <-- FALSE: Recompute standard slack average_converged = checkConvergenceGpu( - iter, d_x_avg_, d_y_avg_, d_ax_avg_, d_aty_avg_, - params_.tolerance, average_results, "[A-GPU]", - d_dSlackPosAvg_, d_dSlackNegAvg_, - false); // <-- FALSE: Recompute standard slack + iter, d_x_avg_, d_y_avg_, d_ax_avg_, d_aty_avg_, params_.tolerance, + average_results, "[A-GPU]", d_dSlackPosAvg_, d_dSlackNegAvg_, + false); // <-- FALSE: Recompute standard slack } #else @@ -842,17 +889,17 @@ bool PDLPSolver::runConvergenceCheckAndRestart(size_t iter, std::vector ATy_pdhg(lp_.num_col_); linalg::Ax(lp_, x_next_, Ax_pdhg); linalg::ATy(lp_, y_next_, ATy_pdhg); - current_converged = checkConvergence(iter, x_next_, y_next_, Ax_pdhg, ATy_pdhg, - params_.tolerance, current_results, "[L]", - dSlackPos_, dSlackNeg_); + current_converged = checkConvergence( + iter, x_next_, y_next_, Ax_pdhg, ATy_pdhg, params_.tolerance, + current_results, "[L]", dSlackPos_, dSlackNeg_); } else { - current_converged = checkConvergence(iter, x_current_, y_current_, Ax_cache_, ATy_cache_, - params_.tolerance, current_results, "[L]", - dSlackPos_, dSlackNeg_); + current_converged = checkConvergence( + iter, x_current_, y_current_, Ax_cache_, ATy_cache_, params_.tolerance, + current_results, "[L]", dSlackPos_, dSlackNeg_); } average_converged = checkConvergence(iter, x_avg_, y_avg_, Ax_avg_, ATy_avg_, - params_.tolerance, average_results, "[A]", - dSlackPosAvg_, dSlackNegAvg_); + params_.tolerance, average_results, + "[A]", dSlackPosAvg_, dSlackNegAvg_); #endif // Determine whether to log iterations @@ -860,47 +907,57 @@ bool PDLPSolver::runConvergenceCheckAndRestart(size_t iter, double time_now = highs_timer_p->read(); iteration_log = time_now > last_logger_time + kHipdlpLoggerFrequency; if (iteration_log) { - logger_.print_iteration_stats(iter, average_results, current_eta_, time_now); + logger_.print_iteration_stats(iter, average_results, current_eta_, + time_now); last_logger_time = time_now; } // 3. Handle Convergence Success if (current_converged || average_converged) { bool prefer_avg = average_converged; - + // Copy result to output vectors #ifdef CUPDLP_GPU double* src_x = prefer_avg ? d_x_avg_ - : (params_.use_halpern_restart ? d_pdhg_primal_ : d_x_current_); + : (params_.use_halpern_restart ? d_pdhg_primal_ + : d_x_current_); double* src_y = prefer_avg ? d_y_avg_ - : (params_.use_halpern_restart ? d_pdhg_dual_ : d_y_current_); + : (params_.use_halpern_restart ? d_pdhg_dual_ + : d_y_current_); double* src_sp = prefer_avg ? d_dSlackPosAvg_ : d_dSlackPos_; double* src_sn = prefer_avg ? d_dSlackNegAvg_ : d_dSlackNeg_; - - CUDA_CHECK(cudaMemcpy(output_x.data(), src_x, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); - CUDA_CHECK(cudaMemcpy(output_y.data(), src_y, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToHost)); - CUDA_CHECK(cudaMemcpy(dSlackPos_.data(), src_sp, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); - CUDA_CHECK(cudaMemcpy(dSlackNeg_.data(), src_sn, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToHost)); + + CUDA_CHECK(cudaMemcpy(output_x.data(), src_x, lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(output_y.data(), src_y, lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(dSlackPos_.data(), src_sp, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); + CUDA_CHECK(cudaMemcpy(dSlackNeg_.data(), src_sn, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToHost)); #else if (prefer_avg) { - output_x = x_avg_; - output_y = y_avg_; - dSlackPos_ = dSlackPosAvg_; - dSlackNeg_ = dSlackNegAvg_; + output_x = x_avg_; + output_y = y_avg_; + dSlackPos_ = dSlackPosAvg_; + dSlackNeg_ = dSlackNegAvg_; } else if (params_.use_halpern_restart) { - output_x = x_next_; - output_y = y_next_; + output_x = x_next_; + output_y = y_next_; } else { - output_x = x_current_; - output_y = y_current_; + output_x = x_current_; + output_y = y_current_; } #endif - + final_iter_count_ = iter; results_ = prefer_avg ? average_results : current_results; - logger_.info((prefer_avg ? "Average" : "Current") + std::string(" solution converged")); - + logger_.info((prefer_avg ? "Average" : "Current") + + std::string(" solution converged")); + status = TerminationStatus::OPTIMAL; - return true; // Stop + return true; // Stop } results_ = current_results; @@ -925,12 +982,12 @@ void PDLPSolver::performPdhgStep() { break; case StepSizeStrategy::MALITSKY_POCK: // Note: Error handling for MP failure simplified for brevity - updateIteratesMalitskyPock(false); + updateIteratesMalitskyPock(false); break; case StepSizeStrategy::PID: - assert(111==999); + assert(111 == 999); } - + #if PDLP_PROFILE hipdlpTimerStop(kHipdlpClockIterateUpdate); #endif @@ -943,7 +1000,7 @@ void PDLPSolver::performHalpernPdhgStep(bool is_major, int k_offset) { double primal_step = stepsize_.primal_step; double dual_step = stepsize_.dual_step; double rho = params_.halpern_gamma; - + int current_k = halpern_iteration_ + k_offset; double w = static_cast(current_k) / (current_k + 1.0); @@ -955,8 +1012,10 @@ void PDLPSolver::performHalpernPdhgStep(bool is_major, int k_offset) { halpern_dual_slack_next_valid_ = true; } for (HighsInt i = 0; i < lp_.num_col_; i++) { - double temp = x_current_[i] - primal_step * (lp_.col_cost_[i] - ATy_cache_[i]); - double temp_proj = linalg::project_box(temp, lp_.col_lower_[i], lp_.col_upper_[i]); + double temp = + x_current_[i] - primal_step * (lp_.col_cost_[i] - ATy_cache_[i]); + double temp_proj = + linalg::project_box(temp, lp_.col_lower_[i], lp_.col_upper_[i]); if (is_major) { x_next_[i] = temp_proj; // pdhg_primal (for convergence check) halpern_dual_slack_next_[i] = (temp_proj - temp) / primal_step; @@ -1012,12 +1071,11 @@ void PDLPSolver::accumulateAverages(size_t iter) { // If Halpern, we average the 'current' blended iterate launchKernelUpdateAverages_wrapper(d_x_sum_, d_y_sum_, d_x_current_, d_y_current_, dMeanStepSize, - lp_.num_col_, lp_.num_row_, - gpu_stream_); + lp_.num_col_, lp_.num_row_, gpu_stream_); sum_weights_gpu_ += dMeanStepSize; } else { // If standard, we average the 'next' iterate computed by PDHG - updateAverageIteratesGpu(inner_iter); + updateAverageIteratesGpu(inner_iter); } #else if (params_.use_halpern_restart) { @@ -1043,10 +1101,18 @@ void PDLPSolver::prepareNextIteration() { } else { // Standard PDHG: x_current becomes x_next #ifdef CUPDLP_GPU - CUDA_CHECK(cudaMemcpy(d_x_current_, d_x_next_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_y_current_, d_y_next_, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_ax_current_, d_ax_next_, lp_.num_row_ * sizeof(double), cudaMemcpyDeviceToDevice)); - CUDA_CHECK(cudaMemcpy(d_aty_current_, d_aty_next_, lp_.num_col_ * sizeof(double), cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_x_current_, d_x_next_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_y_current_, d_y_next_, + lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_ax_current_, d_ax_next_, + lp_.num_row_ * sizeof(double), + cudaMemcpyDeviceToDevice)); + CUDA_CHECK(cudaMemcpy(d_aty_current_, d_aty_next_, + lp_.num_col_ * sizeof(double), + cudaMemcpyDeviceToDevice)); #else x_current_ = x_next_; y_current_ = y_next_; @@ -1062,15 +1128,14 @@ void PDLPSolver::performHalpernPdhgStepGpu(bool is_major, int k_offset) { if (is_major) { launchKernelHalpernPrimalMajor_wrapper( d_x_current_, d_pdhg_primal_, d_x_next_ /*reflected*/, - d_aty_current_ /*dual_product*/, d_col_cost_, - d_col_lower_, d_col_upper_, d_primal_step_size_, a_num_cols_, - d_dSlackPos_ /*dual_slack*/, gpu_stream_); + d_aty_current_ /*dual_product*/, d_col_cost_, d_col_lower_, + d_col_upper_, d_primal_step_size_, a_num_cols_, + d_dSlackPos_ /*dual_slack*/, gpu_stream_); } else { launchKernelHalpernPrimalMinor_wrapper( - d_x_current_, d_x_next_ /*reflected*/, - d_aty_current_ /*dual_product*/, d_col_cost_, - d_col_lower_, d_col_upper_, d_primal_step_size_, a_num_cols_, - gpu_stream_); + d_x_current_, d_x_next_ /*reflected*/, d_aty_current_ /*dual_product*/, + d_col_cost_, d_col_lower_, d_col_upper_, d_primal_step_size_, + a_num_cols_, gpu_stream_); } linalgGpuAx(d_x_next_, d_ax_next_); @@ -1079,22 +1144,22 @@ void PDLPSolver::performHalpernPdhgStepGpu(bool is_major, int k_offset) { launchKernelHalpernDualMajor_wrapper( d_y_current_, d_pdhg_dual_, d_y_next_ /*reflected*/, d_ax_next_ /*primal_product*/, d_row_lower_, d_is_equality_row_, - d_dual_step_size_, a_num_rows_, gpu_stream_); + d_dual_step_size_, a_num_rows_, gpu_stream_); } else { launchKernelHalpernDualMinor_wrapper( - d_y_current_, d_y_next_ /*reflected*/, - d_ax_next_ /*primal_product*/, d_row_lower_, d_is_equality_row_, - d_dual_step_size_, a_num_rows_, gpu_stream_); + d_y_current_, d_y_next_ /*reflected*/, d_ax_next_ /*primal_product*/, + d_row_lower_, d_is_equality_row_, d_dual_step_size_, a_num_rows_, + gpu_stream_); } double rho = params_.halpern_gamma; - launchKernelHalpernBlend_wrapper( - d_x_current_, d_x_next_ /*reflected*/, d_x_anchor_, - d_halpern_iteration_, k_offset, rho, a_num_cols_, gpu_stream_); - launchKernelHalpernBlend_wrapper( - d_y_current_, d_y_next_ /*reflected*/, d_y_anchor_, - d_halpern_iteration_, k_offset, rho, a_num_rows_, gpu_stream_); + launchKernelHalpernBlend_wrapper(d_x_current_, d_x_next_ /*reflected*/, + d_x_anchor_, d_halpern_iteration_, k_offset, + rho, a_num_cols_, gpu_stream_); + launchKernelHalpernBlend_wrapper(d_y_current_, d_y_next_ /*reflected*/, + d_y_anchor_, d_halpern_iteration_, k_offset, + rho, a_num_rows_, gpu_stream_); } #endif @@ -1135,12 +1200,12 @@ void PDLPSolver::initialize() { ATy_cache_.resize(lp_.num_col_, 0.0); Ax_next_.resize(lp_.num_row_, 0.0); ATy_next_.resize(lp_.num_col_, 0.0); - + // Member variables added back in HPP: Ax_avg_.resize(lp_.num_row_, 0.0); ATy_avg_.resize(lp_.num_col_, 0.0); K_times_x_diff_.resize(lp_.num_row_, 0.0); - + dSlackPos_.resize(lp_.num_col_, 0.0); dSlackNeg_.resize(lp_.num_col_, 0.0); dSlackPosAvg_.resize(lp_.num_col_, 0.0); @@ -1305,12 +1370,13 @@ void PDLPSolver::computeDualSlacks(const std::vector& dualResidual, if (params_.use_halpern_restart) { // For Halpern current iterate checks, use major-step dual slack: // dual_slack = (proj(temp) - temp) / primal_step. - // For other cases (e.g. average iterate), fall back to sign/bounds projection. - const bool use_cached_halpern_slack = - (&dSlackPos == &dSlackPos_) && - (&dSlackNeg == &dSlackNeg_) && - halpern_dual_slack_next_valid_ && - (halpern_dual_slack_next_.size() == static_cast(lp_.num_col_)); + // For other cases (e.g. average iterate), fall back to sign/bounds + // projection. + const bool use_cached_halpern_slack = (&dSlackPos == &dSlackPos_) && + (&dSlackNeg == &dSlackNeg_) && + halpern_dual_slack_next_valid_ && + (halpern_dual_slack_next_.size() == + static_cast(lp_.num_col_)); double dual_slack = 0.0; if (use_cached_halpern_slack) { @@ -1397,12 +1463,10 @@ PDLPSolver::computeDualityGap(const std::vector& x, for (HighsInt i = 0; i < lp_.num_col_; ++i) { if (lp_.col_lower_[i] > -kHighsInf) { - lT_lambda_plus += - lp_.col_lower_[i] * std::max(0.0, lambda[i]); + lT_lambda_plus += lp_.col_lower_[i] * std::max(0.0, lambda[i]); } if (lp_.col_upper_[i] < kHighsInf) { - uT_lambda_minus += - lp_.col_upper_[i] * std::max(0.0, -lambda[i]); + uT_lambda_minus += lp_.col_upper_[i] * std::max(0.0, -lambda[i]); } } @@ -1414,7 +1478,8 @@ PDLPSolver::computeDualityGap(const std::vector& x, } double duality_gap = abs(dual_objective - cTx); - return std::make_tuple(duality_gap, qTy, lT_lambda_plus, uT_lambda_minus, cTx); + return std::make_tuple(duality_gap, qTy, lT_lambda_plus, uT_lambda_minus, + cTx); } double PDLPSolver::computeDualObjective(const std::vector& y, @@ -1445,9 +1510,10 @@ double PDLPSolver::computeDualObjective(const std::vector& y, } bool PDLPSolver::checkConvergence( - const HighsInt iter, const std::vector& x, const std::vector& y, - const std::vector& ax_vector, const std::vector& aty_vector, - double epsilon, SolverResults& results, const char* type, + const HighsInt iter, const std::vector& x, + const std::vector& y, const std::vector& ax_vector, + const std::vector& aty_vector, double epsilon, + SolverResults& results, const char* type, // Add slack vectors as non-const references std::vector& dSlackPos, std::vector& dSlackNeg) { // computeDualSlacks is now called inside computeDualFeasibility @@ -1526,9 +1592,9 @@ double PDLPSolver::PowerMethod() { // kYanyuPowerMethodDev; kCuPdlpAATPowerMethod; highsLogDev(params_.log_options_, HighsLogType::kInfo, "Power method: %s\n", - power_method == kYanyuPowerMethod ? "Yanyu" - : power_method == kYanyuPowerMethodDev ? "Yanyu dev" - : "CuPdlp-C"); + power_method == kYanyuPowerMethod ? "Yanyu" + : power_method == kYanyuPowerMethodDev ? "Yanyu dev" + : "CuPdlp-C"); // Dev version of Yanyu power method (based on A'A) has // // * First iterate as vector of ones @@ -1646,15 +1712,14 @@ void PDLPSolver::setup(const HighsOptions& options, HighsTimer& timer) { logger_.passHighsLogOptions(options.log_options); logger_.print_header(); highs_timer_p = &timer; - highsLogUser( - options.log_options, HighsLogType::kInfo, - "Using HiPDLP first order PDLP solver on a %s\n", + highsLogUser(options.log_options, HighsLogType::kInfo, + "Using HiPDLP first order PDLP solver on a %s\n", #ifdef CUPDLP_GPU - "GPU" + "GPU" #else - "CPU: performance may be disappointing!" + "CPU: performance may be disappointing!" #endif - ); + ); #ifdef CUPDLP_GPU HighsInt n_devices = 0; cudaGetDeviceCount(&n_devices); @@ -1685,7 +1750,7 @@ void PDLPSolver::setup(const HighsOptions& options, HighsTimer& timer) { // params.eta = 0; Not set in parse_options_file // params.omega = 0; Not set in parse_options_file params_.tolerance = options.pdlp_optimality_tolerance; - if (options.kkt_tolerance != kDefaultKktTolerance) + if (options.kkt_tolerance != kDefaultKktTolerance) params_.tolerance = options.kkt_tolerance; params_.max_iterations = options.pdlp_iteration_limit; params_.device_type = Device::CPU; @@ -1728,8 +1793,7 @@ void PDLPSolver::setup(const HighsOptions& options, HighsTimer& timer) { } else if (options.pdlp_step_size_strategy == kPdlpStepSizeStrategyMalitskyPock) { params_.step_size_strategy = StepSizeStrategy::MALITSKY_POCK; - } else if (options.pdlp_step_size_strategy == - kPdlpStepSizeStrategyPid) { + } else if (options.pdlp_step_size_strategy == kPdlpStepSizeStrategyPid) { params_.step_size_strategy = StepSizeStrategy::PID; } } @@ -1830,7 +1894,7 @@ void PDLPSolver::initializeStepSizes() { if (params_.step_size_strategy == StepSizeStrategy::FIXED || params_.step_size_strategy == StepSizeStrategy::PID) { // Use power method for fixed step size - //const double op_norm_sq = PowerMethod(); + // const double op_norm_sq = PowerMethod(); double op_norm_sq = PowerMethod(); stepsize_.power_method_lambda = op_norm_sq; @@ -1842,10 +1906,12 @@ void PDLPSolver::initializeStepSizes() { stepsize_.primal_step = base_step / params_.omega; stepsize_.dual_step = base_step * params_.omega; - highsLogDev(params_.log_options_, HighsLogType::kInfo, - "Initial step sizes from power method lambda = %g: primal step= " - "%g; dual step = %g, eta = %g, omega = %g\n", - op_norm_sq, stepsize_.primal_step, stepsize_.dual_step, params_.eta, params_.omega); + highsLogDev( + params_.log_options_, HighsLogType::kInfo, + "Initial step sizes from power method lambda = %g: primal step= " + "%g; dual step = %g, eta = %g, omega = %g\n", + op_norm_sq, stepsize_.primal_step, stepsize_.dual_step, params_.eta, + params_.omega); } else { // Use matrix infinity norm for adaptive step size // Compute infinity norm of matrix elements @@ -1866,9 +1932,9 @@ void PDLPSolver::initializeStepSizes() { stepsize_.dual_step = stepsize_.primal_step * stepsize_.beta; highsLogDev(params_.log_options_, HighsLogType::kInfo, - "Initial step sizes from matrix inf-norm = %g: primal = %g; " - "dual = %g\n", - mat_elem_norm_inf, stepsize_.primal_step, stepsize_.dual_step); + "Initial step sizes from matrix inf-norm = %g: primal = %g; " + "dual = %g\n", + mat_elem_norm_inf, stepsize_.primal_step, stepsize_.dual_step); } } @@ -1876,10 +1942,10 @@ void PDLPSolver::updatePrimalWeightAtRestart(const SolverResults& results) { // Compute delta = pdhg_iterate - anchor double primal_dist = 0.0; double dual_dist = 0.0; - + #ifdef CUPDLP_GPU primal_dist = computeDiffNormCuBLAS(d_pdhg_primal_, d_x_anchor_, a_num_cols_); - dual_dist = computeDiffNormCuBLAS(d_pdhg_dual_, d_y_anchor_, a_num_rows_); + dual_dist = computeDiffNormCuBLAS(d_pdhg_dual_, d_y_anchor_, a_num_rows_); #else // CPU version for (HighsInt i = 0; i < lp_.num_col_; ++i) { @@ -1895,26 +1961,24 @@ void PDLPSolver::updatePrimalWeightAtRestart(const SolverResults& results) { #endif double rel_primal = results.primal_feasibility / (1.0 + unscaled_rhs_norm_); - double rel_dual = results.dual_feasibility / (1.0 + unscaled_c_norm_); + double rel_dual = results.dual_feasibility / (1.0 + unscaled_c_norm_); double ratio_infeas = (rel_primal > 0.0) ? (rel_dual / rel_primal) : 1e300; const bool silent = true; // cuPDLPx-style weight update (PID control) - if (primal_dist > 1e-16 && dual_dist > 1e-16 && - primal_dist < 1e12 && dual_dist < 1e12 && - ratio_infeas > 1e-8 && ratio_infeas < 1e8) { - - double error = std::log(dual_dist) - std::log(primal_dist) - std::log(primal_weight_); - - primal_weight_error_sum_ = params_.i_smooth * primal_weight_error_sum_ + error; + if (primal_dist > 1e-16 && dual_dist > 1e-16 && primal_dist < 1e12 && + dual_dist < 1e12 && ratio_infeas > 1e-8 && ratio_infeas < 1e8) { + double error = + std::log(dual_dist) - std::log(primal_dist) - std::log(primal_weight_); + + primal_weight_error_sum_ = + params_.i_smooth * primal_weight_error_sum_ + error; double delta_error = error - primal_weight_last_error_; - - primal_weight_ *= std::exp( - params_.k_p * error + - params_.k_i * primal_weight_error_sum_ + - params_.k_d * delta_error - ); - + + primal_weight_ *= + std::exp(params_.k_p * error + params_.k_i * primal_weight_error_sum_ + + params_.k_d * delta_error); + primal_weight_last_error_ = error; if (!silent) logger_.info("Primal weight updated: " + std::to_string(primal_weight_)); @@ -1924,21 +1988,25 @@ void PDLPSolver::updatePrimalWeightAtRestart(const SolverResults& results) { primal_weight_error_sum_ = 0.0; primal_weight_last_error_ = 0.0; if (!silent) - logger_.info("Weight update failed (bad norms/ratio), reverted to best: " + - std::to_string(primal_weight_)); + logger_.info( + "Weight update failed (bad norms/ratio), reverted to best: " + + std::to_string(primal_weight_)); } // Track best weight - double gap = (rel_primal > 0.0 && rel_dual > 0.0) ? std::abs(std::log10(rel_dual / rel_primal)) : best_primal_dual_residual_gap_; + double gap = (rel_primal > 0.0 && rel_dual > 0.0) + ? std::abs(std::log10(rel_dual / rel_primal)) + : best_primal_dual_residual_gap_; if (gap < best_primal_dual_residual_gap_) { best_primal_dual_residual_gap_ = gap; best_primal_weight_ = primal_weight_; } - double eta = std::sqrt(stepsize_.primal_step * stepsize_.dual_step); // invariant + double eta = + std::sqrt(stepsize_.primal_step * stepsize_.dual_step); // invariant stepsize_.beta = primal_weight_ * primal_weight_; stepsize_.primal_step = eta / primal_weight_; - stepsize_.dual_step = eta * primal_weight_; + stepsize_.dual_step = eta * primal_weight_; params_.omega = primal_weight_; restart_scheme_.UpdateBeta(stepsize_.beta); } @@ -2071,16 +2139,14 @@ void PDLPSolver::updateIteratesAdaptive() { d_x_current_, // Input (Base) d_aty_current_, // Input (Base ATy) d_col_cost_, d_col_lower_, d_col_upper_, - primal_step_update, a_num_cols_, - gpu_stream_); + primal_step_update, a_num_cols_, gpu_stream_); linalgGpuAx(d_x_next_, d_ax_next_); launchKernelUpdateY_wrapper(d_y_next_, // Output (Trial) d_y_current_, // Input (Base) d_ax_current_, // Input (Base Ax) d_ax_next_, // Input (Trial Ax) d_row_lower_, d_is_equality_row_, - dual_step_update, a_num_rows_, - gpu_stream_); + dual_step_update, a_num_rows_, gpu_stream_); linalgGpuATy(d_y_next_, d_aty_next_); double movement = computeMovementGpu(d_x_next_, d_x_current_, d_y_next_, d_y_current_); @@ -2316,8 +2382,8 @@ void PDLPSolver::setupGpu() { CUDA_CHECK(cudaMemcpy(d_a_row_ptr_, h_a_row_ptr.data(), (a_num_rows_ + 1) * sizeof(HighsInt), cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMemcpy(d_a_col_ind_, h_a_col_ind.data(), a_nnz_ * sizeof(HighsInt), - cudaMemcpyHostToDevice)); + CUDA_CHECK(cudaMemcpy(d_a_col_ind_, h_a_col_ind.data(), + a_nnz_ * sizeof(HighsInt), cudaMemcpyHostToDevice)); CUDA_CHECK(cudaMemcpy(d_a_val_, h_a_val.data(), a_nnz_ * sizeof(double), cudaMemcpyHostToDevice)); @@ -2343,10 +2409,10 @@ void PDLPSolver::setupGpu() { a_nnz_ * sizeof(HighsInt), cudaMemcpyHostToDevice)); CUDA_CHECK(cudaMemcpy(d_at_val_, h_at_val.data(), a_nnz_ * sizeof(double), cudaMemcpyHostToDevice)); - CUDA_CHECK(cudaMalloc(&d_halpern_iteration_, sizeof(int))); - CUDA_CHECK(cudaMemset(d_halpern_iteration_,0, sizeof(int))); + CUDA_CHECK(cudaMalloc(&d_halpern_iteration_, sizeof(int))); + CUDA_CHECK(cudaMemset(d_halpern_iteration_, 0, sizeof(int))); CUDA_CHECK(cudaMalloc(&d_primal_step_size_, sizeof(double))); - CUDA_CHECK(cudaMalloc(&d_dual_step_size_, sizeof(double))); + CUDA_CHECK(cudaMalloc(&d_dual_step_size_, sizeof(double))); // Create AT descriptor with swapped dimensions CUSPARSE_CHECK(cusparseCreateCsr( @@ -2440,7 +2506,7 @@ void PDLPSolver::setupGpu() { &spmv_buffer_size_ax_)); CUDA_CHECK(cudaMalloc(&d_spmv_buffer_ax_, spmv_buffer_size_ax_)); - CUSPARSE_CHECK(cusparseSpMV_preprocess( + CUSPARSE_CHECK(cusparseSpMV_preprocess( cusparse_handle_, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, mat_a_csr_, vec_x_desc_, &beta, vec_ax_desc_, CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG2, d_spmv_buffer_ax_)); @@ -2460,10 +2526,10 @@ void PDLPSolver::setupGpu() { &spmv_buffer_size_aty_)); CUDA_CHECK(cudaMalloc(&d_spmv_buffer_aty_, spmv_buffer_size_aty_)); - CUSPARSE_CHECK(cusparseSpMV_preprocess( - cusparse_handle_, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, - mat_a_T_csr_, vec_y_desc_, &beta, vec_aty_desc_, CUDA_R_64F, - CUSPARSE_SPMV_CSR_ALG2, d_spmv_buffer_aty_)); + CUSPARSE_CHECK(cusparseSpMV_preprocess( + cusparse_handle_, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, mat_a_T_csr_, + vec_y_desc_, &beta, vec_aty_desc_, CUDA_R_64F, CUSPARSE_SPMV_CSR_ALG2, + d_spmv_buffer_aty_)); CUSPARSE_CHECK(cusparseDestroyDnVec(vec_y)); CUSPARSE_CHECK(cusparseDestroyDnVec(vec_aty)); @@ -2502,8 +2568,8 @@ void PDLPSolver::setupGpu() { CUDA_CHECK(cudaMalloc(&d_buffer2_, max_size * sizeof(double))); highsLogDev(params_.log_options_, HighsLogType::kInfo, - "GPU setup complete. Matrix A (CSR) and A^T (CSR) transferred " - "to device.\n"); + "GPU setup complete. Matrix A (CSR) and A^T (CSR) transferred " + "to device.\n"); } void PDLPSolver::cleanupGpu() { @@ -2642,10 +2708,10 @@ bool PDLPSolver::checkConvergenceGpu(const HighsInt iter, const double* d_x, double* d_slackNeg_out, bool use_halpern_slack) { launchCheckConvergenceKernels_wrapper( - d_convergence_results_, d_slackPos_out, d_slackNeg_out, d_x, d_y, d_ax, d_aty, - d_col_cost_, d_row_lower_, d_col_lower_, d_col_upper_, d_is_equality_row_, - d_col_scale_, d_row_scale_, lp_.num_col_, lp_.num_row_, use_halpern_slack, - gpu_stream_); + d_convergence_results_, d_slackPos_out, d_slackNeg_out, d_x, d_y, d_ax, + d_aty, d_col_cost_, d_row_lower_, d_col_lower_, d_col_upper_, + d_is_equality_row_, d_col_scale_, d_row_scale_, lp_.num_col_, + lp_.num_row_, use_halpern_slack, gpu_stream_); // copy 4 doubles back to CPU diff --git a/highs/pdlp/hipdlp/pdhg.hpp b/highs/pdlp/hipdlp/pdhg.hpp index 939c997e26..59f7a99509 100644 --- a/highs/pdlp/hipdlp/pdhg.hpp +++ b/highs/pdlp/hipdlp/pdhg.hpp @@ -21,8 +21,8 @@ #include #include -#include #include +#include #include "linalg.hpp" #include "logger.hpp" @@ -83,10 +83,10 @@ class PDLPSolver { void passLp(const HighsLp* lp) { original_lp_ = lp; } void preprocessLp(); void scaleProblem(); - + // Main entry point void solve(std::vector& x, std::vector& y); - + void unscaleSolution(std::vector& x, std::vector& y); PostSolveRetcode postprocess(HighsSolution& solution); void logSummary(); @@ -110,23 +110,22 @@ class PDLPSolver { // --- Core Algorithm Logic --- void initialize(); void solveReturn(const TerminationStatus term_code); - + // Returns true if solver should terminate, false if it should continue // Updates 'status' if returning true. - bool runConvergenceCheckAndRestart(size_t iter, - std::vector& output_x, + bool runConvergenceCheckAndRestart(size_t iter, std::vector& output_x, std::vector& output_y, TerminationStatus& status); void performPdhgStep(); void performHalpernStep(); void accumulateAverages(size_t iter); - void prepareNextIteration(); // Swaps pointers/vectors + void prepareNextIteration(); // Swaps pointers/vectors // --- Convergence & Math Helpers --- double PowerMethod(); void initializeStepSizes(); - + // Convergence checks bool checkConvergence(const HighsInt iter, const std::vector& x, const std::vector& y, @@ -136,53 +135,71 @@ class PDLPSolver { std::vector& dSlackPos, std::vector& dSlackNeg); - void computeAverageIterate(std::vector& ax_avg, std::vector& aty_avg); - + void computeAverageIterate(std::vector& ax_avg, + std::vector& aty_avg); + // Step updates - std::vector updateX(const std::vector& x, const std::vector& aty, double primal_step); - std::vector updateY(const std::vector& y, const std::vector& ax, const std::vector& ax_next, double dual_step); + std::vector updateX(const std::vector& x, + const std::vector& aty, + double primal_step); + std::vector updateY(const std::vector& y, + const std::vector& ax, + const std::vector& ax_next, + double dual_step); void updateIteratesFixed(); void updateIteratesAdaptive(); bool updateIteratesMalitskyPock(bool first_malitsky_iteration); // Restart helpers void computeStepSizeRatio(PrimalDualParams& working_params); - void applyHalpernAveraging(std::vector& x, std::vector& y, std::vector& ax, std::vector& aty); - void updateAverageIterates(const std::vector& x, const std::vector& y, HighsInt inner_iter); + void applyHalpernAveraging(std::vector& x, std::vector& y, + std::vector& ax, std::vector& aty); + void updateAverageIterates(const std::vector& x, + const std::vector& y, HighsInt inner_iter); void performHalpernPdhgStep(bool is_major, int k_offset); void updatePrimalWeightAtRestart(const SolverResults& results); // Feasibility calculations double computePrimalFeasibility(const std::vector& Ax_vector); - double computeDualFeasibility(const std::vector& ATy_vector, std::vector& dSlackPos, std::vector& dSlackNeg); - void computeDualSlacks(const std::vector& dualResidual, std::vector& dSlackPos, std::vector& dSlackNeg); - double computeDualObjective(const std::vector& y, const std::vector& dSlackPos, const std::vector& dSlackNeg); - std::vector computeLambda(const std::vector& y, const std::vector& ATy_vector); + double computeDualFeasibility(const std::vector& ATy_vector, + std::vector& dSlackPos, + std::vector& dSlackNeg); + void computeDualSlacks(const std::vector& dualResidual, + std::vector& dSlackPos, + std::vector& dSlackNeg); + double computeDualObjective(const std::vector& y, + const std::vector& dSlackPos, + const std::vector& dSlackNeg); + std::vector computeLambda(const std::vector& y, + const std::vector& ATy_vector); std::tuple computeDualityGap( const std::vector& x, const std::vector& y, const std::vector& lambda); // Other utilities void printConstraintInfo(); - bool CheckNumericalStability(const std::vector& delta_x, const std::vector& delta_y); - double computeMovement(const std::vector& delta_primal, const std::vector& delta_dual); - double computeNonlinearity(const std::vector& delta_primal, const std::vector& delta_aty); + bool CheckNumericalStability(const std::vector& delta_x, + const std::vector& delta_y); + double computeMovement(const std::vector& delta_primal, + const std::vector& delta_dual); + double computeNonlinearity(const std::vector& delta_primal, + const std::vector& delta_aty); double computeFixedPointError(); // --- Data Members --- HighsLp lp_; const HighsLp* original_lp_ = nullptr; HighsLp unscaled_processed_lp_; - + PrimalDualParams params_; StepSizeConfig stepsize_; Logger logger_; HighsLogOptions log_options_; SolverResults results_; - + Scaling scaling_; RestartScheme restart_scheme_; - + // Problem dimensions & metadata HighsInt original_num_col_ = 0; HighsInt num_eq_rows_ = 0; @@ -199,21 +216,21 @@ class PDLPSolver { std::vector x_avg_, y_avg_; std::vector x_sum_, y_sum_; std::vector x_at_last_restart_, y_at_last_restart_; - std::vector reflected_x_, reflected_y_; // For over-relaxed Halpern - std::vector x_anchor_, y_anchor_; // For Halpern + std::vector reflected_x_, reflected_y_; // For over-relaxed Halpern + std::vector x_anchor_, y_anchor_; // For Halpern // Caches std::vector Ax_cache_, ATy_cache_; std::vector Ax_next_, ATy_next_; // Average Caches - std::vector Ax_avg_, ATy_avg_; + std::vector Ax_avg_, ATy_avg_; // Residual Caches std::vector K_times_x_diff_; std::vector dSlackPos_, dSlackNeg_; std::vector dSlackPosAvg_, dSlackNegAvg_; std::vector halpern_dual_slack_next_; bool halpern_dual_slack_next_valid_ = false; - double initial_fpe_ = 0.0; // For Halpern restart + double initial_fpe_ = 0.0; // For Halpern restart double fpe_ = 0.0; // Scalars @@ -229,7 +246,8 @@ class PDLPSolver { double best_primal_weight_ = 1.0; double primal_weight_error_sum_ = 0.0; double primal_weight_last_error_ = 0.0; - double best_primal_dual_residual_gap_ = std::numeric_limits::infinity(); + double best_primal_dual_residual_gap_ = + std::numeric_limits::infinity(); HighsTimer* highs_timer_p; double last_logger_time = -kHighsInf; @@ -251,16 +269,16 @@ class PDLPSolver { cudaStream_t gpu_stream_ = nullptr; cusparseHandle_t cusparse_handle_ = nullptr; cublasHandle_t cublas_handle_ = nullptr; - + // Matrix descriptors cusparseSpMatDescr_t mat_a_csr_ = nullptr; cusparseSpMatDescr_t mat_a_T_csr_ = nullptr; - + // Device Pointers HighsInt *d_a_row_ptr_ = nullptr, *d_a_col_ind_ = nullptr; - double *d_a_val_ = nullptr; + double* d_a_val_ = nullptr; HighsInt *d_at_row_ptr_ = nullptr, *d_at_col_ind_ = nullptr; - double *d_at_val_ = nullptr; + double* d_at_val_ = nullptr; int* d_halpern_iteration_ = nullptr; double* d_primal_step_size_ = nullptr; double* d_dual_step_size_ = nullptr; @@ -273,29 +291,34 @@ class PDLPSolver { double *d_ax_current_ = nullptr, *d_aty_current_ = nullptr; double *d_ax_next_ = nullptr, *d_aty_next_ = nullptr; double *d_ax_avg_ = nullptr, *d_aty_avg_ = nullptr; - + double *d_x_at_last_restart_ = nullptr, *d_y_at_last_restart_ = nullptr; double *d_x_anchor_ = nullptr, *d_y_anchor_ = nullptr; double *d_pdhg_primal_ = nullptr, *d_pdhg_dual_ = nullptr; double* d_delta_x_ = nullptr; double* d_delta_y_ = nullptr; double* d_AT_delta_y_ = nullptr; - - double *d_col_cost_ = nullptr, *d_col_lower_ = nullptr, *d_col_upper_ = nullptr; - double *d_row_lower_ = nullptr, *d_col_scale_ = nullptr, *d_row_scale_ = nullptr; - bool *d_is_equality_row_ = nullptr; + + double *d_col_cost_ = nullptr, *d_col_lower_ = nullptr, + *d_col_upper_ = nullptr; + double *d_row_lower_ = nullptr, *d_col_scale_ = nullptr, + *d_row_scale_ = nullptr; + bool* d_is_equality_row_ = nullptr; // Scratch / Output buffers - double *d_convergence_results_ = nullptr; + double* d_convergence_results_ = nullptr; double *d_dSlackPos_ = nullptr, *d_dSlackNeg_ = nullptr; double *d_dSlackPosAvg_ = nullptr, *d_dSlackNegAvg_ = nullptr; - double *d_buffer_ = nullptr, *d_buffer2_ = nullptr; // General purpose - double *d_x_temp_diff_norm_result_ = nullptr, *d_y_temp_diff_norm_result_ = nullptr; + double *d_buffer_ = nullptr, *d_buffer2_ = nullptr; // General purpose + double *d_x_temp_diff_norm_result_ = nullptr, + *d_y_temp_diff_norm_result_ = nullptr; // SpMV buffers - void *d_spmv_buffer_ax_ = nullptr; size_t spmv_buffer_size_ax_ = 0; - void *d_spmv_buffer_aty_ = nullptr; size_t spmv_buffer_size_aty_ = 0; - + void* d_spmv_buffer_ax_ = nullptr; + size_t spmv_buffer_size_ax_ = 0; + void* d_spmv_buffer_aty_ = nullptr; + size_t spmv_buffer_size_aty_ = 0; + // Vector Descriptors cusparseDnVecDescr_t vec_x_desc_ = nullptr, vec_y_desc_ = nullptr; cusparseDnVecDescr_t vec_ax_desc_ = nullptr, vec_aty_desc_ = nullptr; @@ -305,16 +328,20 @@ class PDLPSolver { void cleanupGpu(); void linalgGpuAx(const double* d_x_in, double* d_ax_out); void linalgGpuATy(const double* d_y_in, double* d_aty_out); - bool checkConvergenceGpu(const HighsInt iter, const double* d_x, const double* d_y, - const double* d_ax, const double* d_aty, - double epsilon, SolverResults& results, - const char* type, double* d_slackPos_out, double* d_slackNeg_out, + bool checkConvergenceGpu(const HighsInt iter, const double* d_x, + const double* d_y, const double* d_ax, + const double* d_aty, double epsilon, + SolverResults& results, const char* type, + double* d_slackPos_out, double* d_slackNeg_out, bool use_halpern_slack); void computeStepSizeRatioGpu(PrimalDualParams& working_params); void updateAverageIteratesGpu(HighsInt inner_iter); void computeAverageIterateGpu(); - double computeMovementGpu(const double* d_x_new, const double* d_x_old, const double* d_y_new, const double* d_y_old); - double computeNonlinearityGpu(const double* d_x_new, const double* d_x_old, const double* d_aty_new, const double* d_aty_old); + double computeMovementGpu(const double* d_x_new, const double* d_x_old, + const double* d_y_new, const double* d_y_old); + double computeNonlinearityGpu(const double* d_x_new, const double* d_x_old, + const double* d_aty_new, + const double* d_aty_old); double computeDiffNormCuBLAS(const double* d_a, const double* d_b, int n); void performHalpernPdhgStepGpu(bool is_major, int k_offset); double computeFixedPointErrorGpu(); @@ -323,8 +350,9 @@ class PDLPSolver { void launchKernelUpdateX(double primal_step); void launchKernelUpdateY(double dual_step); void launchKernelUpdateAverages(double weight); - void launchKernelScaleVector(double* d_out, const double* d_in, double scale, HighsInt n); + void launchKernelScaleVector(double* d_out, const double* d_in, double scale, + HighsInt n); #endif }; -#endif // PDHG_HPP +#endif // PDHG_HPP diff --git a/highs/pdlp/hipdlp/scaling.cc b/highs/pdlp/hipdlp/scaling.cc index 28eea0dd90..4cc678695c 100644 --- a/highs/pdlp/hipdlp/scaling.cc +++ b/highs/pdlp/hipdlp/scaling.cc @@ -79,21 +79,21 @@ void Scaling::scaleProblem() { if (params_->use_ruiz_scaling) { highsLogDev(params_->log_options_, HighsLogType::kInfo, - "Applying Ruiz scaling...\n"); + "Applying Ruiz scaling...\n"); ApplyRuizScaling(); is_scaled_ = true; } if (params_->use_pc_scaling) { highsLogDev(params_->log_options_, HighsLogType::kInfo, - "Applying Pock-Chambolle scaling...\n"); + "Applying Pock-Chambolle scaling...\n"); ApplyPockChambolleScaling(); is_scaled_ = true; } if (params_->use_l2_scaling) { highsLogDev(params_->log_options_, HighsLogType::kInfo, - "Applying L2 norm scaling...\n"); + "Applying L2 norm scaling...\n"); ApplyL2Scaling(); is_scaled_ = true; } @@ -132,7 +132,7 @@ void Scaling::ApplyRuizScaling() { highsLogUser(params_->log_options_, HighsLogType::kError, "Currently only support infinity norm for Ruiz scaling\n"); assert(ruiz_norm_ok); - } + } if (params_->ruiz_norm == INFINITY) { // For infinity norm, find max absolute value in each row for (HighsInt col = 0; col < lp_->num_col_; ++col) { @@ -168,9 +168,7 @@ void Scaling::ApplyRuizScaling() { } void Scaling::ApplyPockChambolleScaling() { - const bool pc_alpha_ok = - 0.0 <= params_->pc_alpha && - params_->pc_alpha <= 2.0; + const bool pc_alpha_ok = 0.0 <= params_->pc_alpha && params_->pc_alpha <= 2.0; if (!pc_alpha_ok) { highsLogUser(params_->log_options_, HighsLogType::kError, "PC alpha must be in [0, 2]\n"); diff --git a/highs/pdlp/hipdlp/scaling.hpp b/highs/pdlp/hipdlp/scaling.hpp index 7c2dbe29b3..414701b86a 100644 --- a/highs/pdlp/hipdlp/scaling.hpp +++ b/highs/pdlp/hipdlp/scaling.hpp @@ -35,7 +35,7 @@ class Scaling { debug_pdlp_data_ = debug_pdlp_data; }; // Not used - // void LogMatrixNorms(const std::string& stage); + // void LogMatrixNorms(const std::string& stage); // // Get scaling vectors (for unscaling solution later) bool IsScaled() const { return is_scaled_; } @@ -65,7 +65,8 @@ class Scaling { const std::vector& row_scaling); // Compute norm of a vector based on norm type - double ComputeNorm(const double* values, HighsInt size, double norm_type) const; + double ComputeNorm(const double* values, HighsInt size, + double norm_type) const; FILE* debug_pdlp_log_file_; DebugPdlpData* debug_pdlp_data_; }; diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 8edfa6e518..7fae95b87a 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -5662,13 +5662,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 @@ -5678,6 +5681,7 @@ HPresolve::Result HPresolve::initialRowAndColPresolve( HPRESOLVE_CHECKED_CALL(colPresolve(postsolve_stack, col)); changedColFlag[col] = false; } + analysis_.presolveTimerStop(kPresolveClockInitialCol); return checkLimits(postsolve_stack); } @@ -5687,15 +5691,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); @@ -5775,7 +5789,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, @@ -5809,7 +5825,9 @@ HPresolve::Result HPresolve::presolve(HighsPostsolveStack& postsolve_stack) { report(); } + analysis_.presolveTimerStart(kPresolveClockFastLoop); HPRESOLVE_CHECKED_CALL(fastPresolveLoop(postsolve_stack)); + analysis_.presolveTimerStop(kPresolveClockFastLoop); storeCurrentProblemSize(); @@ -5822,14 +5840,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(numNz))); @@ -5837,12 +5860,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; } @@ -5851,7 +5871,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_); @@ -5859,12 +5883,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; @@ -5877,7 +5905,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(); @@ -5910,7 +5940,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_); @@ -5919,11 +5951,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; } @@ -6410,6 +6446,9 @@ HPresolve::Result HPresolve::removeDependentEquations( matrix.index_.reserve(maxCapacity); std::vector eqSet(matrix.num_col_); + std::vector row_count; + row_count.assign(model->num_col_, 0); + HighsInt i = 0; for (const std::pair& p : equations) { HighsInt eq = p.second; @@ -6417,8 +6456,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,6 +6470,9 @@ HPresolve::Result HPresolve::removeDependentEquations( matrix.start_[i] = matrix.value_.size(); } + HighsInt num_nonzero_columns = 0; + for (HighsInt iCol = 0; iCol < model->num_col_; iCol++) + if (row_count[iCol]) num_nonzero_columns++; std::vector colSet(matrix.num_col_); std::iota(colSet.begin(), colSet.end(), 0); HFactor factor; @@ -6464,7 +6508,7 @@ HPresolve::Result HPresolve::removeDependentEquations( time_taken); analysis_.logging_on_ = logging_on; if (logging_on) - analysis_.stopPresolveRuleLog(kPresolveRuleDependentFreeCols); + analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations); return Result::kOk; } else { double pct_off_timeout = @@ -6494,12 +6538,28 @@ HPresolve::Result HPresolve::removeDependentEquations( num_fictitious_rows_skipped++; } } - if (!silent) + 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 (limit = %.2fs)\n", + static_cast(matrix.num_col_), matrix.num_col_ == 1 ? "" : "s", + static_cast(num_nonzero_columns), + static_cast(model->num_col_), num_nonzero_columns == 1 ? "" : "s", + static_cast(matrix.numNz()), matrix.numNz() == 1 ? "" : "s", + static_cast(num_removed_row), num_removed_row == 1 ? "" : "s", + static_cast(num_removed_nz), num_removed_nz == 1 ? "" : "s", + time_taken, time_limit); highsLogUser(options->log_options, HighsLogType::kInfo, - "Dependent equations search removed %d rows and %d nonzeros " - "in %.2fs (limit = %.2fs)\n", + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", + model->model_name_.c_str(), static_cast(matrix.num_col_), + static_cast(num_nonzero_columns), + static_cast(model->num_col_), + static_cast(matrix.numNz()), static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken, time_limit); + static_cast(num_removed_nz), time_taken); + } if (num_fictitious_rows_skipped) highsLogDev(options->log_options, HighsLogType::kInfo, ", avoiding %d fictitious rows", diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index a96a4b0d5e..1cd1897613 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -13,15 +13,15 @@ void HPresolveAnalysis::setup(const HighsLp* model_, const HighsOptions* options_, const HighsInt& numDeletedRows_, const HighsInt& numDeletedCols_, - const bool silent, - HighsTimer* timer) { + const bool silent, HighsTimer* timer) { model = model_; options = options_; numDeletedRows = &numDeletedRows_; numDeletedCols = &numDeletedCols_; timer_ = timer; - analyse_presolve_time_ = kHighsAnalysisLevelPresolveTime & options->highs_analysis_level; + analyse_presolve_time_ = + kHighsAnalysisLevelPresolveTime & options->highs_analysis_level; if (analyse_presolve_time_) { HighsTimerClock clock; clock.timer_pointer_ = timer_; @@ -250,7 +250,8 @@ bool HPresolveAnalysis::analysePresolveRuleLog(const bool report) { return true; } -void HPresolveAnalysis::presolveTimerStart(const HighsInt presolve_clock) const { +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); @@ -265,6 +266,5 @@ void HPresolveAnalysis::presolveTimerStop(const HighsInt presolve_clock) const { void HPresolveAnalysis::reportPresolveTimer() { if (!analyse_presolve_time_) return; PresolveTimer presolve_timer; - presolve_timer.reportPresolveCoreClock(presolve_clocks_); + presolve_timer.reportPresolveCoreClock(model->model_name_, presolve_clocks_); } - diff --git a/highs/presolve/HPresolveAnalysis.h b/highs/presolve/HPresolveAnalysis.h index b6ba39fb00..e61eec96ca 100644 --- a/highs/presolve/HPresolveAnalysis.h +++ b/highs/presolve/HPresolveAnalysis.h @@ -15,10 +15,8 @@ class HPresolveAnalysis { public: - HPresolveAnalysis() - : timer_(nullptr), - analyse_presolve_time_(false) {} - + HPresolveAnalysis() : timer_(nullptr), analyse_presolve_time_(false) {} + HighsTimer* timer_; const HighsLp* model; const HighsOptions* options; @@ -40,7 +38,7 @@ class HPresolveAnalysis { HighsInt num_deleted_rows0_; HighsInt num_deleted_cols0_; HighsPresolveLog presolve_log_; - + HighsTimerClock presolve_clocks_; bool analyse_presolve_time_; diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h index 23622b54a4..8d33223158 100644 --- a/highs/presolve/PresolveTimer.h +++ b/highs/presolve/PresolveTimer.h @@ -15,6 +15,23 @@ enum iClockPresolve { kPresolveClockTotal = 0, kPresolveClockPresolve, + kPresolveClockInitial, + kPresolveClockInitialRow, + kPresolveClockInitialChangeColBounds, + kPresolveClockInitialCol, + kPresolveClockFastLoop, + kPresolveClockFastLoopRowSingletons, + kPresolveClockFastLoopColSingletons, + kPresolveClockFastLoopDoubletonEquations, + kPresolveClockFastLoopChangedRows, + kPresolveClockFastLoopChangedCols, + kPresolveClockAggregator, + kPresolveClockSparsify, + kPresolveClockParallelRowsAndCols, + kPresolveClockDependentEquations, + kPresolveClockDependentFreeCol, + kPresolveClockShrinkProblem, + // kPresolveClock@, kNumPresolveClock //!< Number of PRESOLVE clocks }; @@ -29,13 +46,40 @@ class PresolveTimer { 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[kPresolveClockInitialChangeColBounds] = + timer_pointer->clock_def("Initial change col bounds"); + 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 presolve_clock_list, - const HighsTimerClock& presolve_timer_clock, - const HighsInt kPresolveClockIdeal = kPresolveClockTotal, - const double tolerance_percent_report_ = -1) { + bool reportPresolveClockList( + const char* grepStamp, const std::vector 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& clock = presolve_timer_clock.clock_; @@ -53,11 +97,66 @@ class PresolveTimer { grepStamp, clockList, ideal_sum_time, tolerance_percent_report); }; - void reportPresolveCoreClock(const HighsTimerClock& presolve_timer_clock) { + void csvPresolveClockList(const std::string& grep_query, + const std::string& model_name, + const std::vector 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& 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 presolve_clock_list{ - kPresolveClockPresolve}; - reportPresolveClockList("PresolveCore_", presolve_clock_list, presolve_timer_clock, - kPresolveClockTotal); + // kPresolveClockInitial, + kPresolveClockInitialRow, kPresolveClockInitialChangeColBounds, + 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); }; }; From 4518f70d9df68ebd4b97439053f0fa45a363a4b2 Mon Sep 17 00:00:00 2001 From: JAJHall Date: Wed, 25 Mar 2026 23:48:25 +0000 Subject: [PATCH 05/16] Aded heuristic to perform dependent equation presolve rule; removed presolve-time unit test as it's not silent --- check/TestPresolve.cpp | 10 ------ highs/presolve/HPresolve.cpp | 49 +++++++++++++++++++--------- highs/presolve/HPresolveAnalysis.cpp | 46 +++++++++++++++----------- highs/presolve/PresolveTimer.h | 6 +--- 4 files changed, 61 insertions(+), 50 deletions(-) diff --git a/check/TestPresolve.cpp b/check/TestPresolve.cpp index f8478a0fe2..8af6618f34 100644 --- a/check/TestPresolve.cpp +++ b/check/TestPresolve.cpp @@ -910,13 +910,3 @@ TEST_CASE("presolve-issue-2874", "[highs_test_presolve]") { REQUIRE(highs.presolve() == HighsStatus::kOk); REQUIRE(highs.getModelPresolveStatus() == HighsPresolveStatus::kInfeasible); } - -TEST_CASE("presolve-time", "[highs_test_presolve]") { - std::string model_file = - std::string(HIGHS_DIR) + "/check/instances/shell.mps"; - Highs highs; - // highs.setOptionValue("output_flag", dev_run); - highs.readModel(model_file); - highs.setOptionValue("highs_analysis_level", 256); - REQUIRE(highs.presolve() == HighsStatus::kOk); -} diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 7fae95b87a..da039dacd7 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6433,11 +6433,19 @@ 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_[0] = 0; @@ -6473,6 +6481,17 @@ HPresolve::Result HPresolve::removeDependentEquations( HighsInt num_nonzero_columns = 0; for (HighsInt iCol = 0; iCol < model->num_col_; iCol++) if (row_count[iCol]) num_nonzero_columns++; + HighsInt num_variables = num_nonzero_columns; + // + // Although two duplicated equations in any number of variables lead + // to dependency, it's not worth looking for dependent equations if + // there are few relative to the number of variables + // + // The following heuristic picks up all Mittelmann problems with + // meaningful reductions + if (num_equations < 0.5 * num_variables) return returnOk(); + // + // Identify any dependent equations std::vector colSet(matrix.num_col_); std::iota(colSet.begin(), colSet.end(), 0); HFactor factor; @@ -6539,11 +6558,19 @@ HPresolve::Result HPresolve::removeDependentEquations( } } if (!silent) { + highsLogUser(options->log_options, HighsLogType::kInfo, + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", + model->model_name_.c_str(), static_cast(matrix.num_col_), + static_cast(num_nonzero_columns), + static_cast(model->num_col_), + static_cast(matrix.numNz()), + static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); 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 (limit = %.2fs)\n", + "in %.2fs (limit = %.2fs)", static_cast(matrix.num_col_), matrix.num_col_ == 1 ? "" : "s", static_cast(num_nonzero_columns), static_cast(model->num_col_), num_nonzero_columns == 1 ? "" : "s", @@ -6551,25 +6578,15 @@ HPresolve::Result HPresolve::removeDependentEquations( static_cast(num_removed_row), num_removed_row == 1 ? "" : "s", static_cast(num_removed_nz), num_removed_nz == 1 ? "" : "s", time_taken, time_limit); - highsLogUser(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(matrix.num_col_), - static_cast(num_nonzero_columns), - static_cast(model->num_col_), - static_cast(matrix.numNz()), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); } if (num_fictitious_rows_skipped) highsLogDev(options->log_options, HighsLogType::kInfo, - ", avoiding %d fictitious rows", - static_cast(num_fictitious_rows_skipped)); + ", avoiding %d fictitious row%s", + static_cast(num_fictitious_rows_skipped), + num_fictitious_rows_skipped == 1 ? "" : "s"); highsLogDev(options->log_options, HighsLogType::kInfo, "\n"); - analysis_.logging_on_ = logging_on; - if (logging_on) - analysis_.stopPresolveRuleLog(kPresolveRuleDependentEquations); - return Result::kOk; + return returnOk(); } HPresolve::Result HPresolve::removeDependentFreeCols( diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index 1cd1897613..79cb39d430 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -21,7 +21,8 @@ void HPresolveAnalysis::setup(const HighsLp* model_, timer_ = timer; analyse_presolve_time_ = - kHighsAnalysisLevelPresolveTime & options->highs_analysis_level; + kHighsAnalysisLevelPresolveTime & options->highs_analysis_level && + !model_->isMip(); if (analyse_presolve_time_) { HighsTimerClock clock; clock.timer_pointer_ = timer_; @@ -32,21 +33,30 @@ void HPresolveAnalysis::setup(const HighsLp* model_, 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 = pow(2, 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++) { @@ -56,13 +66,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 diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h index 8d33223158..138c6fb418 100644 --- a/highs/presolve/PresolveTimer.h +++ b/highs/presolve/PresolveTimer.h @@ -17,7 +17,6 @@ enum iClockPresolve { kPresolveClockPresolve, kPresolveClockInitial, kPresolveClockInitialRow, - kPresolveClockInitialChangeColBounds, kPresolveClockInitialCol, kPresolveClockFastLoop, kPresolveClockFastLoopRowSingletons, @@ -48,8 +47,6 @@ class PresolveTimer { clock[kPresolveClockPresolve] = timer_pointer->clock_def("Presolve"); clock[kPresolveClockInitial] = timer_pointer->clock_def("Initial"); clock[kPresolveClockInitialRow] = timer_pointer->clock_def("Initial row"); - clock[kPresolveClockInitialChangeColBounds] = - timer_pointer->clock_def("Initial change col bounds"); clock[kPresolveClockInitialCol] = timer_pointer->clock_def("Initial col"); clock[kPresolveClockFastLoop] = timer_pointer->clock_def("Fast loop"); clock[kPresolveClockFastLoopRowSingletons] = @@ -137,8 +134,7 @@ class PresolveTimer { const HighsTimerClock& presolve_timer_clock) { const std::vector presolve_clock_list{ // kPresolveClockInitial, - kPresolveClockInitialRow, kPresolveClockInitialChangeColBounds, - kPresolveClockInitialCol, + kPresolveClockInitialRow, kPresolveClockInitialCol, // kPresolveClockFastLoop, kPresolveClockFastLoopRowSingletons, kPresolveClockFastLoopColSingletons, From 8bc59ca2d3ba24363e72b85c80f20e2422fa55cd Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Thu, 26 Mar 2026 13:02:16 +0000 Subject: [PATCH 06/16] Use of std::pow no unambiguous? --- highs/presolve/HPresolve.cpp | 13 +++++++------ highs/presolve/HPresolveAnalysis.cpp | 2 +- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index da039dacd7..b6bb8d7a8d 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -1740,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(kPresolveRuleProbing)))); return Result::kStopped; } diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index 79cb39d430..748814c66a 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -38,7 +38,7 @@ void HPresolveAnalysis::setup(const HighsLp* model_, highsLogUser(options->log_options, HighsLogType::kInfo, "Permitted suppression of presolve rules via " "presolve_rule_off option:\n"); - HighsInt bit = pow(2, kPresolveRuleFirstAllowOff); + HighsInt bit = pow(int(2), int(kPresolveRuleFirstAllowOff)); for (HighsInt rule_type = kPresolveRuleFirstAllowOff; rule_type < kPresolveRuleCount; rule_type++) { // This is a rule that can be switched off From ff138f5334322646f587ba3c9af886457efa5ffb Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Thu, 26 Mar 2026 14:58:07 +0000 Subject: [PATCH 07/16] Reduced kMaxDependentEquationsTime to 10 --- highs/io/HighsIO.cpp | 5 ++ highs/io/HighsIO.h | 1 + highs/presolve/HPresolve.cpp | 90 ++++++++++++++++++---------- highs/presolve/HPresolveAnalysis.cpp | 3 +- highs/presolve/PresolveTimer.h | 2 + 5 files changed, 67 insertions(+), 34 deletions(-) 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/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index b6bb8d7a8d..2e949966c5 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6448,13 +6448,13 @@ HPresolve::Result HPresolve::removeDependentEquations( 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 eqSet(matrix.num_col_); + std::vector eqSet(num_equations); std::vector row_count; row_count.assign(model->num_col_, 0); @@ -6479,10 +6479,12 @@ HPresolve::Result HPresolve::removeDependentEquations( matrix.start_[i] = matrix.value_.size(); } - HighsInt num_nonzero_columns = 0; + // 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_nonzero_columns++; - HighsInt num_variables = num_nonzero_columns; + if (row_count[iCol]) num_variables++; + HighsInt num_nz = matrix.numNz(); // // Although two duplicated equations in any number of variables lead // to dependency, it's not worth looking for dependent equations if @@ -6490,10 +6492,21 @@ HPresolve::Result HPresolve::removeDependentEquations( // // The following heuristic picks up all Mittelmann problems with // meaningful reductions - if (num_equations < 0.5 * num_variables) return returnOk(); + const bool not_consider_dependent_equations = + num_equations < 0.5 * num_variables; + const bool silent = silentLog(); + if (!silent) + highsLogUser(options->log_options, HighsLogType::kInfo, + "%sonsidering dependency of %d equation%s in %d variable%s " + "with %d nonzero%s\n", + not_consider_dependent_equations ? "Not c" : "C", + int(num_equations), highsIntToPlural(num_equations).c_str(), + int(num_variables), highsIntToPlural(num_variables).c_str(), + int(num_nz), highsIntToPlural(num_nz).c_str()); + if (not_consider_dependent_equations) return returnOk(); // // Identify any dependent equations - std::vector colSet(matrix.num_col_); + std::vector colSet(num_equations); std::iota(colSet.begin(), colSet.end(), 0); HFactor factor; factor.setup(matrix, colSet); @@ -6506,30 +6519,42 @@ 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 = 10; + 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(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) { + /* + highsLogUser(options->log_options, HighsLogType::kInfo, + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g,Terminated\n", + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), + static_cast(model->num_col_), + static_cast(num_nz), + static_cast(num_removed_row), + static_cast(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(kPresolveRuleDependentEquations); - return Result::kOk; + } + return returnOk(); } else { double pct_off_timeout = 1e2 * std::fabs(time_taken - time_limit) / time_limit; @@ -6543,10 +6568,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]]; @@ -6559,32 +6580,35 @@ HPresolve::Result HPresolve::removeDependentEquations( } } if (!silent) { + /* highsLogUser(options->log_options, HighsLogType::kInfo, "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(matrix.num_col_), - static_cast(num_nonzero_columns), + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), static_cast(model->num_col_), - static_cast(matrix.numNz()), + static_cast(num_nz), static_cast(num_removed_row), static_cast(num_removed_nz), time_taken); + */ 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 (limit = %.2fs)", - static_cast(matrix.num_col_), matrix.num_col_ == 1 ? "" : "s", - static_cast(num_nonzero_columns), - static_cast(model->num_col_), num_nonzero_columns == 1 ? "" : "s", - static_cast(matrix.numNz()), matrix.numNz() == 1 ? "" : "s", - static_cast(num_removed_row), num_removed_row == 1 ? "" : "s", - static_cast(num_removed_nz), num_removed_nz == 1 ? "" : "s", - time_taken, time_limit); + static_cast(num_equations), + highsIntToPlural(num_equations).c_str(), + static_cast(num_variables), static_cast(model->num_col_), + highsIntToPlural(num_variables).c_str(), static_cast(num_nz), + num_nz == 1 ? "" : "s", static_cast(num_removed_row), + highsIntToPlural(num_removed_row).c_str(), + static_cast(num_removed_nz), + highsIntToPlural(num_removed_nz).c_str(), time_taken, time_limit); } if (num_fictitious_rows_skipped) highsLogDev(options->log_options, HighsLogType::kInfo, ", avoiding %d fictitious row%s", static_cast(num_fictitious_rows_skipped), - num_fictitious_rows_skipped == 1 ? "" : "s"); + highsIntToPlural(num_fictitious_rows_skipped).c_str()); highsLogDev(options->log_options, HighsLogType::kInfo, "\n"); return returnOk(); diff --git a/highs/presolve/HPresolveAnalysis.cpp b/highs/presolve/HPresolveAnalysis.cpp index 748814c66a..8ca3362847 100644 --- a/highs/presolve/HPresolveAnalysis.cpp +++ b/highs/presolve/HPresolveAnalysis.cpp @@ -38,7 +38,8 @@ void HPresolveAnalysis::setup(const HighsLp* model_, highsLogUser(options->log_options, HighsLogType::kInfo, "Permitted suppression of presolve rules via " "presolve_rule_off option:\n"); - HighsInt bit = pow(int(2), int(kPresolveRuleFirstAllowOff)); + HighsInt bit = + std::pow(int(2), static_cast(kPresolveRuleFirstAllowOff)); for (HighsInt rule_type = kPresolveRuleFirstAllowOff; rule_type < kPresolveRuleCount; rule_type++) { // This is a rule that can be switched off diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h index 138c6fb418..a8b992c43c 100644 --- a/highs/presolve/PresolveTimer.h +++ b/highs/presolve/PresolveTimer.h @@ -147,12 +147,14 @@ class PresolveTimer { }; 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); + */ }; }; From ea86ff40d9e07e4c7f114aa0cb0f1848c32e5c8b Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Thu, 26 Mar 2026 15:02:35 +0000 Subject: [PATCH 08/16] kMaxDependentEquationsTime=10 is too small: try kMaxDependentEquationsTime=100 --- highs/presolve/HPresolve.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 2e949966c5..76250476ba 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6519,7 +6519,7 @@ 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 kMaxDependentEquationsTime = 10; + const double kMaxDependentEquationsTime = 100; const double time_limit = std::max( 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); From 447c66be14b9ac11a699040d6d83b1feeaa81dc8 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Thu, 26 Mar 2026 16:17:33 +0000 Subject: [PATCH 09/16] kMaxDependentEquationsTime back to 1000; formatted --- highs/presolve/HPresolve.cpp | 38 +++++++++++++++------------------- highs/presolve/PresolveTimer.h | 2 -- 2 files changed, 17 insertions(+), 23 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 76250476ba..9d3fc03a28 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6519,7 +6519,7 @@ 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 kMaxDependentEquationsTime = 100; + const double kMaxDependentEquationsTime = 1000; const double time_limit = std::max( 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); @@ -6539,16 +6539,14 @@ HPresolve::Result HPresolve::removeDependentEquations( if (build_return == kBuildKernelReturnTimeout) { // HFactor::build has timed out, so just return if (!silent) { - /* - highsLogUser(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g,Terminated\n", - model->model_name_.c_str(), static_cast(num_equations), - static_cast(num_variables), - static_cast(model->num_col_), - static_cast(num_nz), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); - */ + 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(num_equations), + static_cast(num_variables), static_cast(model->num_col_), + static_cast(num_nz), static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); highsLogUser(options->log_options, HighsLogType::kInfo, "Dependent equations search terminated after %.3gs due to " "expected time exceeding limit\n", @@ -6580,16 +6578,14 @@ HPresolve::Result HPresolve::removeDependentEquations( } } if (!silent) { - /* - highsLogUser(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(num_equations), - static_cast(num_variables), - static_cast(model->num_col_), - static_cast(num_nz), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); - */ + if (options->log_dev_level > 0) + highsLogUser(options->log_options, HighsLogType::kInfo, + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), + static_cast(model->num_col_), static_cast(num_nz), + static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); highsLogUser( options->log_options, HighsLogType::kInfo, "Search of %d equation%s with %d / %d variable%s and %d nonzero%s " diff --git a/highs/presolve/PresolveTimer.h b/highs/presolve/PresolveTimer.h index a8b992c43c..138c6fb418 100644 --- a/highs/presolve/PresolveTimer.h +++ b/highs/presolve/PresolveTimer.h @@ -147,14 +147,12 @@ class PresolveTimer { }; 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); - */ }; }; From 3e58ddeedfd46247f79d0ca4007ce6a5310c0f6d Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Fri, 27 Mar 2026 11:06:27 +0000 Subject: [PATCH 10/16] Now to correct kernel_num_el --- highs/presolve/HPresolve.cpp | 4 +- highs/util/HFactor.cpp | 76 ++++++++++++++++++++++++++++++++++-- highs/util/HFactor.h | 3 +- 3 files changed, 76 insertions(+), 7 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 9d3fc03a28..5197e97073 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6503,7 +6503,7 @@ HPresolve::Result HPresolve::removeDependentEquations( int(num_equations), highsIntToPlural(num_equations).c_str(), int(num_variables), highsIntToPlural(num_variables).c_str(), int(num_nz), highsIntToPlural(num_nz).c_str()); - if (not_consider_dependent_equations) return returnOk(); + // if (not_consider_dependent_equations) return returnOk(); // // Identify any dependent equations std::vector colSet(num_equations); @@ -6519,7 +6519,7 @@ 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 kMaxDependentEquationsTime = 1000; + const double kMaxDependentEquationsTime = 100; const double time_limit = std::max( 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index 31a16af254..e7b67db546 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; @@ -821,14 +825,17 @@ void HFactor::buildSimple() { HighsInt mr_countX = 0; // Determine the number of entries in the kernel kernel_num_el = 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; + check_nwork++; } } mr_index.resize(mr_countX); @@ -841,6 +848,7 @@ void HFactor::buildSimple() { mc_count_a.assign(mc_dim, 0); mc_count_n.assign(mc_dim, 0); HighsInt MCcountX = 0; + num_active_nz_ = 0; for (HighsInt i = 0; i < nwork; i++) { HighsInt iCol = iwork[i]; mc_var[iCol] = b_var[iCol]; @@ -853,9 +861,12 @@ 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); + num_active_nz_++; } else { + // Above the active part of the kernel colStoreN(iCol, iRow, value); } } @@ -865,6 +876,10 @@ void HFactor::buildSimple() { build_synthetic_tick += (num_row + nwork + MCcountX) * 40 + mr_countX * 20; // Record the kernel dimension kernel_dim = nwork; + if (num_row == num_basic) { + assert(kernel_dim == check_nwork); + assert(kernel_num_el == num_active_nz_ + kernel_dim); + } assert((HighsInt)this->refactor_info_.pivot_row.size() == num_basic - nwork); } @@ -884,12 +899,37 @@ HighsInt HFactor::buildKernel() { const bool check_for_timeout = this->time_limit_ < kHighsInf; HighsInt search_k = 0; + auto numActiveNz = [&]() { + // Determine the number of indices in the active submatrix + HighsInt num_active_nz = 0; + const HighsInt max_count = max(num_row, num_basic); + for (HighsInt count = 2; count <= max_count; count++) { + // Column count cannot exceed the number of rows + if (count <= num_row) { + for (HighsInt j = col_link_first[count]; j != -1; j = col_link_next[j]) + num_active_nz += mc_count_a[j]; + } + } + return num_active_nz; + }; + + HighsInt num_active_nz = numActiveNz(); + + printf("HFactor::buildKernel Initial num_active_nz = %d\n", int(num_active_nz)); const HighsInt check_nwork = -11; while (nwork-- > 0) { // printf("\nnwork = %d\n", (int)nwork); if (nwork == check_nwork) { reportAsm(); } + HighsInt check_num_active_nz = numActiveNz(); + /* + if (num_active_nz != check_num_active_nz) { + printf("nwork = %d: search_k = %d, num_active_nz = %d != %d = check_num_active_nz\n", + int(nwork), int(search_k), int(num_active_nz), int(check_num_active_nz)); + exit(1); + } + */ // Determine whether to return due to exceeding the time limit if (check_for_timeout && search_k % timer_frequency == 0) { double current_time = build_timer_->read(); @@ -899,14 +939,35 @@ HighsInt HFactor::buildKernel() { average_iteration_time = 0.9 * average_iteration_time + 0.1 * iteration_time; - if (time_difference > this->time_limit_ / 1e3) + if (search_k > 0 && time_difference > this->time_limit_ / 1e3) { + HighsInt prev_timer_frequency = timer_frequency; timer_frequency = std::max(HighsInt(1), timer_frequency / 10); + if (timer_frequency != prev_timer_frequency) { + printf("Timer frequency change from %d to %d since time_difference = %g > %g = time_limit_ / 1e3\n", + int(prev_timer_frequency), int(timer_frequency), + time_difference, this->time_limit_ / 1e3); + } + } + check_num_active_nz = numActiveNz(); + assert(num_active_nz == check_num_active_nz); + HighsInt iterations_left = kernel_dim - search_k + 1; double remaining_time_bound = average_iteration_time * iterations_left; double total_time_bound = current_time + remaining_time_bound; + const bool stop = current_time > this->time_limit_ || + total_time_bound > this->time_limit_; + printf("HFactor::buildKernel search_k = %7d; kernel_dim = %7d; num_active_nz = %7d;" + " Time (Cu = %6.2f; Iter = %6.4f; AvIter = %.6f; Bd = %.6f Lim = %6.2f) Fq = %4d - %s\n", + int(search_k), int(kernel_dim), int(num_active_nz), + current_time, iteration_time, average_iteration_time, + total_time_bound, this->time_limit_, + int(timer_frequency), + stop ? "Stop" : ""); + /* if (current_time > this->time_limit_ || total_time_bound > this->time_limit_) return kBuildKernelReturnTimeout; + */ } /** @@ -1044,7 +1105,8 @@ HighsInt HFactor::buildKernel() { assert(jColPivot < 0); assert(!foundPivot); rank_deficiency = nwork + 1; - highsLogDev(log_options, HighsLogType::kWarning, + // highsLogDev(log_options, HighsLogType::kWarning, + printf( "Factorization identifies rank deficiency of %d\n", (int)rank_deficiency); return rank_deficiency; @@ -1061,6 +1123,7 @@ HighsInt HFactor::buildKernel() { // // Remove the pivot row index from the pivotal column of the // col-wise matrix. Also decreases the column count + 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 @@ -1122,6 +1185,7 @@ HighsInt HFactor::buildKernel() { l_value.push_back(value); mr_count_before[iRow] = mr_count[iRow]; rowDelete(jColPivot, (int)iRow); + num_active_nz--; } l_start.push_back(l_index.size()); fake_fill += 2 * mc_count_a[jColPivot]; @@ -1149,6 +1213,7 @@ HighsInt HFactor::buildKernel() { const HighsInt my_end = my_start + my_count - 1; double my_pivot = colDelete(iCol, iRowPivot); colStoreN(iCol, iRowPivot, my_pivot); + num_active_nz--; // 2.4.2. Elimination on the overlapping part HighsInt nFillin = mwz_column_count; @@ -1167,6 +1232,7 @@ HighsInt HFactor::buildKernel() { mc_value[my_k] = value; } } + num_active_nz -= nCancel; fake_eliminate += mwz_column_count; fake_eliminate += nFillin * 2; @@ -1207,8 +1273,10 @@ 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]); + num_active_nz++; + } } // 2.4.4.3 Fill into the row copy diff --git a/highs/util/HFactor.h b/highs/util/HFactor.h index 8e452846f8..7cd7e2d78d 100644 --- a/highs/util/HFactor.h +++ b/highs/util/HFactor.h @@ -336,7 +336,8 @@ class HFactor { HighsInt invert_num_el; HighsInt kernel_dim; HighsInt kernel_num_el; - + HighsInt num_active_nz_; + /** * Data of the factor */ From 6ceeddb7567f5b8c6bcdda4753fd7fcaeefccac9 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Fri, 27 Mar 2026 23:40:46 +0000 Subject: [PATCH 11/16] Move timeout code down --- highs/presolve/HPresolve.cpp | 28 +++--- highs/simplex/HighsSimplexAnalysis.cpp | 4 +- highs/util/HFactor.cpp | 121 ++++++++++++++----------- 3 files changed, 85 insertions(+), 68 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 5197e97073..06e29e4c5e 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6578,14 +6578,13 @@ HPresolve::Result HPresolve::removeDependentEquations( } } if (!silent) { - if (options->log_dev_level > 0) - highsLogUser(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(num_equations), - static_cast(num_variables), - static_cast(model->num_col_), static_cast(num_nz), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); + highsLogDev(options->log_options, HighsLogType::kInfo, + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), + static_cast(model->num_col_), static_cast(num_nz), + static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); highsLogUser( options->log_options, HighsLogType::kInfo, "Search of %d equation%s with %d / %d variable%s and %d nonzero%s " @@ -6599,14 +6598,13 @@ HPresolve::Result HPresolve::removeDependentEquations( highsIntToPlural(num_removed_row).c_str(), static_cast(num_removed_nz), highsIntToPlural(num_removed_nz).c_str(), time_taken, time_limit); + if (num_fictitious_rows_skipped) + highsLogDev(options->log_options, HighsLogType::kInfo, + ", avoiding %d fictitious row%s", + static_cast(num_fictitious_rows_skipped), + highsIntToPlural(num_fictitious_rows_skipped).c_str()); + highsLogUser(options->log_options, HighsLogType::kInfo, "\n"); } - if (num_fictitious_rows_skipped) - highsLogDev(options->log_options, HighsLogType::kInfo, - ", avoiding %d fictitious row%s", - static_cast(num_fictitious_rows_skipped), - highsIntToPlural(num_fictitious_rows_skipped).c_str()); - highsLogDev(options->log_options, HighsLogType::kInfo, "\n"); - return returnOk(); } diff --git a/highs/simplex/HighsSimplexAnalysis.cpp b/highs/simplex/HighsSimplexAnalysis.cpp index 1d24b7b2fb..dcaf2ae545 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; + (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 e7b67db546..0d1400f063 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -824,7 +824,7 @@ void HFactor::buildSimple() { mr_count.assign(num_row, 0); HighsInt mr_countX = 0; // Determine the number of entries in the kernel - kernel_num_el = 0; + num_active_nz_ = 0; HighsInt check_nwork = 0; for (HighsInt iRow = 0; iRow < num_row; iRow++) { HighsInt count = mr_count_before[iRow]; @@ -834,7 +834,7 @@ void HFactor::buildSimple() { mr_space[iRow] = count * 2; mr_countX += count * 2; rlinkAdd(iRow, count); - kernel_num_el += count + 1; + num_active_nz_ += count; check_nwork++; } } @@ -848,7 +848,6 @@ void HFactor::buildSimple() { mc_count_a.assign(mc_dim, 0); mc_count_n.assign(mc_dim, 0); HighsInt MCcountX = 0; - num_active_nz_ = 0; for (HighsInt i = 0; i < nwork; i++) { HighsInt iCol = iwork[i]; mc_var[iCol] = b_var[iCol]; @@ -864,7 +863,6 @@ void HFactor::buildSimple() { // In the active part of the kernel colInsert(iCol, iRow, value); rowInsert(iCol, iRow); - num_active_nz_++; } else { // Above the active part of the kernel colStoreN(iCol, iRow, value); @@ -874,12 +872,9 @@ 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; - if (num_row == num_basic) { - assert(kernel_dim == check_nwork); - assert(kernel_num_el == num_active_nz_ + kernel_dim); - } + kernel_num_el = num_active_nz_; assert((HighsInt)this->refactor_info_.pivot_row.size() == num_basic - nwork); } @@ -892,46 +887,42 @@ HighsInt HFactor::buildKernel() { const bool progress_report = false; // num_basic != num_row; const HighsInt progress_frequency = 10000; + + HighsInt search_k = 0; + HighsInt num_defer_pivot = 0; // Initial timer frequency: may be reduced if iterations get slow - HighsInt timer_frequency = 100; - double previous_iteration_time = 0; + HighsInt timer_frequency = 10; + HighsInt log_frequency = 1000; + double previous_iteration_time = build_timer_->read(); double average_iteration_time = 0; + HighsInt previous_num_active_nz = num_active_nz_; + double average_num_active_nz_change_rate = 0; const bool check_for_timeout = this->time_limit_ < kHighsInf; - HighsInt search_k = 0; - auto numActiveNz = [&]() { + auto checkNumActiveNz = [&]() { // Determine the number of indices in the active submatrix - HighsInt num_active_nz = 0; - const HighsInt max_count = max(num_row, num_basic); - for (HighsInt count = 2; count <= max_count; count++) { - // Column count cannot exceed the number of rows - if (count <= num_row) { - for (HighsInt j = col_link_first[count]; j != -1; j = col_link_next[j]) - num_active_nz += mc_count_a[j]; - } + HighsInt check_num_active_nz = 0; + for (HighsInt count = 1; count <= num_row; count++) { + for (HighsInt j = col_link_first[count]; j != -1; j = col_link_next[j]) + check_num_active_nz += count; + } + if (num_active_nz_ != check_num_active_nz) { + printf("nwork = %d: search_k = %d, num_active_nz_ = %d != %d = check_num_active_nz\n", + int(nwork), int(search_k), int(num_active_nz_), int(check_num_active_nz)); + assert(num_active_nz_ == check_num_active_nz); + num_active_nz_ = check_num_active_nz; } - return num_active_nz; }; - HighsInt num_active_nz = numActiveNz(); - - printf("HFactor::buildKernel Initial num_active_nz = %d\n", int(num_active_nz)); const HighsInt check_nwork = -11; while (nwork-- > 0) { // printf("\nnwork = %d\n", (int)nwork); if (nwork == check_nwork) { reportAsm(); } - HighsInt check_num_active_nz = numActiveNz(); - /* - if (num_active_nz != check_num_active_nz) { - printf("nwork = %d: search_k = %d, num_active_nz = %d != %d = check_num_active_nz\n", - int(nwork), int(search_k), int(num_active_nz), int(check_num_active_nz)); - exit(1); - } - */ // Determine whether to return due to exceeding the time limit - if (check_for_timeout && search_k % timer_frequency == 0) { + if (check_for_timeout && search_k > 0 && search_k % timer_frequency == 0) { + checkNumActiveNz(); double current_time = build_timer_->read(); double time_difference = current_time - previous_iteration_time; previous_iteration_time = current_time; @@ -939,7 +930,9 @@ HighsInt HFactor::buildKernel() { average_iteration_time = 0.9 * average_iteration_time + 0.1 * iteration_time; - if (search_k > 0 && time_difference > this->time_limit_ / 1e3) { + // Determine whether to reduce the timer frequency + /* + if (time_difference > this->time_limit_ / 1e1) { HighsInt prev_timer_frequency = timer_frequency; timer_frequency = std::max(HighsInt(1), timer_frequency / 10); if (timer_frequency != prev_timer_frequency) { @@ -948,21 +941,46 @@ HighsInt HFactor::buildKernel() { time_difference, this->time_limit_ / 1e3); } } - check_num_active_nz = numActiveNz(); - assert(num_active_nz == check_num_active_nz); - - HighsInt iterations_left = kernel_dim - search_k + 1; + */ + double num_active_nz_change_rate = + static_cast(num_active_nz_ - previous_num_active_nz) / + static_cast(timer_frequency); + previous_num_active_nz = num_active_nz_; + average_num_active_nz_change_rate = + 0.9 * average_num_active_nz_change_rate + 0.1 * 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; + HighsInt iterations_left0 = iterations_left; + HighsInt iterations_left1 = 0; + if (average_num_active_nz_change_rate < 0) { + HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; + iterations_left1 = active_nz_iterations_left; + 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; const bool stop = current_time > this->time_limit_ || total_time_bound > this->time_limit_; - printf("HFactor::buildKernel search_k = %7d; kernel_dim = %7d; num_active_nz = %7d;" - " Time (Cu = %6.2f; Iter = %6.4f; AvIter = %.6f; Bd = %.6f Lim = %6.2f) Fq = %4d - %s\n", - int(search_k), int(kernel_dim), int(num_active_nz), - current_time, iteration_time, average_iteration_time, - total_time_bound, this->time_limit_, - int(timer_frequency), - stop ? "Stop" : ""); + if (search_k % log_frequency == 0) { + printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" + " ActiveNzRate(%6d, Av = %6d);" + " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s\n", + int(search_k), int(kernel_dim), + int(num_defer_pivot), + int(iterations_left0), + int(iterations_left1), + int(iterations_left), + int(num_active_nz_), + int(num_active_nz_change_rate), + int(average_num_active_nz_change_rate), + current_time, 1000*iteration_time, 1000*average_iteration_time, + total_time_bound, this->time_limit_, + int(timer_frequency), + stop ? " - Stop" : ""); + } /* if (current_time > this->time_limit_ || total_time_bound > this->time_limit_) @@ -1123,7 +1141,7 @@ HighsInt HFactor::buildKernel() { // // Remove the pivot row index from the pivotal column of the // col-wise matrix. Also decreases the column count - num_active_nz--; + 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 @@ -1141,6 +1159,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, @@ -1185,7 +1204,7 @@ HighsInt HFactor::buildKernel() { l_value.push_back(value); mr_count_before[iRow] = mr_count[iRow]; rowDelete(jColPivot, (int)iRow); - num_active_nz--; + num_active_nz_--; } l_start.push_back(l_index.size()); fake_fill += 2 * mc_count_a[jColPivot]; @@ -1213,7 +1232,7 @@ HighsInt HFactor::buildKernel() { const HighsInt my_end = my_start + my_count - 1; double my_pivot = colDelete(iCol, iRowPivot); colStoreN(iCol, iRowPivot, my_pivot); - num_active_nz--; + num_active_nz_--; // 2.4.2. Elimination on the overlapping part HighsInt nFillin = mwz_column_count; @@ -1232,7 +1251,7 @@ HighsInt HFactor::buildKernel() { mc_value[my_k] = value; } } - num_active_nz -= nCancel; + num_active_nz_ -= nCancel; fake_eliminate += mwz_column_count; fake_eliminate += nFillin * 2; @@ -1275,7 +1294,7 @@ HighsInt HFactor::buildKernel() { HighsInt iRow = mwz_column_index[i]; if (mwz_column_mark[iRow]) { colInsert(iCol, iRow, -my_pivot * mwz_column_array[iRow]); - num_active_nz++; + num_active_nz_++; } } From 8969c4c4babd9e8fd055d642b5c1cba426e28bfb Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Fri, 27 Mar 2026 23:45:39 +0000 Subject: [PATCH 12/16] Move timeout code down --- highs/util/HFactor.cpp | 139 +++++++++++++++++++++-------------------- 1 file changed, 71 insertions(+), 68 deletions(-) diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index 0d1400f063..9bb25b2e5c 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -896,6 +896,7 @@ HighsInt HFactor::buildKernel() { double previous_iteration_time = build_timer_->read(); double average_iteration_time = 0; HighsInt previous_num_active_nz = num_active_nz_; + HighsInt previous_search_k = search_k; double average_num_active_nz_change_rate = 0; const bool check_for_timeout = this->time_limit_ < kHighsInf; @@ -920,74 +921,6 @@ HighsInt HFactor::buildKernel() { if (nwork == check_nwork) { reportAsm(); } - // Determine whether to return due to exceeding the time limit - if (check_for_timeout && search_k > 0 && search_k % timer_frequency == 0) { - checkNumActiveNz(); - 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; - - // Determine whether to reduce the timer frequency - /* - if (time_difference > this->time_limit_ / 1e1) { - HighsInt prev_timer_frequency = timer_frequency; - timer_frequency = std::max(HighsInt(1), timer_frequency / 10); - if (timer_frequency != prev_timer_frequency) { - printf("Timer frequency change from %d to %d since time_difference = %g > %g = time_limit_ / 1e3\n", - int(prev_timer_frequency), int(timer_frequency), - time_difference, this->time_limit_ / 1e3); - } - } - */ - double num_active_nz_change_rate = - static_cast(num_active_nz_ - previous_num_active_nz) / - static_cast(timer_frequency); - previous_num_active_nz = num_active_nz_; - average_num_active_nz_change_rate = - 0.9 * average_num_active_nz_change_rate + 0.1 * 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; - HighsInt iterations_left0 = iterations_left; - HighsInt iterations_left1 = 0; - if (average_num_active_nz_change_rate < 0) { - HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; - iterations_left1 = active_nz_iterations_left; - 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; - const bool stop = current_time > this->time_limit_ || - total_time_bound > this->time_limit_; - if (search_k % log_frequency == 0) { - printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" - " ActiveNzRate(%6d, Av = %6d);" - " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s\n", - int(search_k), int(kernel_dim), - int(num_defer_pivot), - int(iterations_left0), - int(iterations_left1), - int(iterations_left), - int(num_active_nz_), - int(num_active_nz_change_rate), - int(average_num_active_nz_change_rate), - current_time, 1000*iteration_time, 1000*average_iteration_time, - total_time_bound, this->time_limit_, - int(timer_frequency), - stop ? " - Stop" : ""); - } - /* - if (current_time > this->time_limit_ || - total_time_bound > this->time_limit_) - return kBuildKernelReturnTimeout; - */ - } - /** * 1. Search for the pivot */ @@ -1115,6 +1048,76 @@ HighsInt HFactor::buildKernel() { fake_search += count; } } + // Determine whether to return due to exceeding the time limit + if (check_for_timeout && search_k > 0 && (search_k % timer_frequency == 0 ||iRowPivot < 0)) { + HighsInt dl_search_k = search_k - previous_search_k; + previous_search_k = search_k; + checkNumActiveNz(); + 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 * dl_search_k); + average_iteration_time = + 0.9 * average_iteration_time + 0.1 * iteration_time; + + // Determine whether to reduce the timer frequency + /* + if (time_difference > this->time_limit_ / 1e1) { + HighsInt prev_timer_frequency = timer_frequency; + timer_frequency = std::max(HighsInt(1), timer_frequency / 10); + if (timer_frequency != prev_timer_frequency) { + printf("Timer frequency change from %d to %d since time_difference = %g > %g = time_limit_ / 1e3\n", + int(prev_timer_frequency), int(timer_frequency), + time_difference, this->time_limit_ / 1e3); + } + } + */ + double num_active_nz_change_rate = + static_cast(num_active_nz_ - previous_num_active_nz) / + static_cast(dl_search_k); + previous_num_active_nz = num_active_nz_; + average_num_active_nz_change_rate = + 0.9 * average_num_active_nz_change_rate + 0.1 * 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; + HighsInt iterations_left0 = iterations_left; + HighsInt iterations_left1 = 0; + if (average_num_active_nz_change_rate < 0) { + HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; + iterations_left1 = active_nz_iterations_left; + 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; + const bool stop = current_time > this->time_limit_ || + total_time_bound > this->time_limit_; + if (search_k % log_frequency == 0) { + printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" + " ActiveNzRate(%6d, Av = %6d);" + " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s\n", + int(search_k), int(kernel_dim), + int(num_defer_pivot), + int(iterations_left0), + int(iterations_left1), + int(iterations_left), + int(num_active_nz_), + int(num_active_nz_change_rate), + int(average_num_active_nz_change_rate), + current_time, 1000*iteration_time, 1000*average_iteration_time, + total_time_bound, this->time_limit_, + int(timer_frequency), + stop ? " - Stop" : ""); + } + /* + if (current_time > this->time_limit_ || + total_time_bound > this->time_limit_) + return kBuildKernelReturnTimeout; + */ + } + // 1.4. If we found nothing: tell singular if (iRowPivot < 0) { // To detect the absence of a pivot, it should be sufficient From 1516fba75ebc21c346e75e2df6146524fe156f34 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Sat, 28 Mar 2026 18:51:33 +0000 Subject: [PATCH 13/16] Use resizeIntDataShift --- highs/presolve/HPresolve.cpp | 2 +- highs/util/HFactor.cpp | 173 +++++++++++++++++------------------ 2 files changed, 86 insertions(+), 89 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 06e29e4c5e..7f32427744 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6519,7 +6519,7 @@ 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 kMaxDependentEquationsTime = 100; + const double kMaxDependentEquationsTime = 1000; const double time_limit = std::max( 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index 9bb25b2e5c..2e4b67f789 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -891,36 +891,92 @@ HighsInt HFactor::buildKernel() { HighsInt search_k = 0; HighsInt num_defer_pivot = 0; // Initial timer frequency: may be reduced if iterations get slow + const HighsInt max_timer_frequency = 1000; HighsInt timer_frequency = 10; HighsInt log_frequency = 1000; double previous_iteration_time = build_timer_->read(); double average_iteration_time = 0; HighsInt previous_num_active_nz = num_active_nz_; - HighsInt previous_search_k = search_k; double average_num_active_nz_change_rate = 0; + bool resize_data_shift = false; const bool check_for_timeout = this->time_limit_ < kHighsInf; - auto checkNumActiveNz = [&]() { - // Determine the number of indices in the active submatrix - HighsInt check_num_active_nz = 0; - for (HighsInt count = 1; count <= num_row; count++) { - for (HighsInt j = col_link_first[count]; j != -1; j = col_link_next[j]) - check_num_active_nz += count; - } - if (num_active_nz_ != check_num_active_nz) { - printf("nwork = %d: search_k = %d, num_active_nz_ = %d != %d = check_num_active_nz\n", - int(nwork), int(search_k), int(num_active_nz_), int(check_num_active_nz)); - assert(num_active_nz_ == check_num_active_nz); - num_active_nz_ = check_num_active_nz; - } + auto resizeIntDataShift = [&](std::vector 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; }; - + 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 > 0 && search_k % timer_frequency == 0) { + 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); + const double mu0 = (1.0 * timer_frequency) / (10.0 * max_timer_frequency); + assert(mu0 <= 0.2); + const double mu1 = 1.0 - mu0; + if (!resize_data_shift) { + average_iteration_time = + mu1 * average_iteration_time + mu0 * iteration_time; + + // Determine whether to reduce the timer frequency + if (time_difference > this->time_limit_ / 1e1) + timer_frequency = std::max(HighsInt(1), timer_frequency / 10); + } + double num_active_nz_change_rate = + static_cast(num_active_nz_ - previous_num_active_nz) / + static_cast(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; + HighsInt iterations_left0 = iterations_left; + HighsInt iterations_left1 = 0; + if (average_num_active_nz_change_rate < 0) { + HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; + iterations_left1 = active_nz_iterations_left; + 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; + const bool stop = current_time > this->time_limit_ || + total_time_bound > this->time_limit_; + if (search_k % log_frequency == 0 || stop || resize_data_shift) { + printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" + " ActiveNzRate(%6d, Av = %6d);" + " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s%s\n", + int(search_k), int(kernel_dim), + int(num_defer_pivot), + int(iterations_left0), + int(iterations_left1), + int(iterations_left), + int(num_active_nz_), + int(num_active_nz_change_rate), + int(average_num_active_nz_change_rate), + current_time, 1000*iteration_time, 1000*average_iteration_time, + total_time_bound, this->time_limit_, + int(timer_frequency), + resize_data_shift ? " - ReZ-shift " : "", + stop ? " - Stop" : ""); + } + /* + if (current_time > this->time_limit_ || + total_time_bound > this->time_limit_) + return kBuildKernelReturnTimeout; + */ + resize_data_shift = false; + } /** * 1. Search for the pivot */ @@ -1048,76 +1104,6 @@ HighsInt HFactor::buildKernel() { fake_search += count; } } - // Determine whether to return due to exceeding the time limit - if (check_for_timeout && search_k > 0 && (search_k % timer_frequency == 0 ||iRowPivot < 0)) { - HighsInt dl_search_k = search_k - previous_search_k; - previous_search_k = search_k; - checkNumActiveNz(); - 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 * dl_search_k); - average_iteration_time = - 0.9 * average_iteration_time + 0.1 * iteration_time; - - // Determine whether to reduce the timer frequency - /* - if (time_difference > this->time_limit_ / 1e1) { - HighsInt prev_timer_frequency = timer_frequency; - timer_frequency = std::max(HighsInt(1), timer_frequency / 10); - if (timer_frequency != prev_timer_frequency) { - printf("Timer frequency change from %d to %d since time_difference = %g > %g = time_limit_ / 1e3\n", - int(prev_timer_frequency), int(timer_frequency), - time_difference, this->time_limit_ / 1e3); - } - } - */ - double num_active_nz_change_rate = - static_cast(num_active_nz_ - previous_num_active_nz) / - static_cast(dl_search_k); - previous_num_active_nz = num_active_nz_; - average_num_active_nz_change_rate = - 0.9 * average_num_active_nz_change_rate + 0.1 * 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; - HighsInt iterations_left0 = iterations_left; - HighsInt iterations_left1 = 0; - if (average_num_active_nz_change_rate < 0) { - HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; - iterations_left1 = active_nz_iterations_left; - 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; - const bool stop = current_time > this->time_limit_ || - total_time_bound > this->time_limit_; - if (search_k % log_frequency == 0) { - printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" - " ActiveNzRate(%6d, Av = %6d);" - " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s\n", - int(search_k), int(kernel_dim), - int(num_defer_pivot), - int(iterations_left0), - int(iterations_left1), - int(iterations_left), - int(num_active_nz_), - int(num_active_nz_change_rate), - int(average_num_active_nz_change_rate), - current_time, 1000*iteration_time, 1000*average_iteration_time, - total_time_bound, this->time_limit_, - int(timer_frequency), - stop ? " - Stop" : ""); - } - /* - if (current_time > this->time_limit_ || - total_time_bound > this->time_limit_) - return kBuildKernelReturnTimeout; - */ - } - // 1.4. If we found nothing: tell singular if (iRowPivot < 0) { // To detect the absence of a pivot, it should be sufficient @@ -1126,13 +1112,11 @@ HighsInt HFactor::buildKernel() { assert(jColPivot < 0); assert(!foundPivot); rank_deficiency = nwork + 1; - // highsLogDev(log_options, HighsLogType::kWarning, - printf( + highsLogDev(log_options, HighsLogType::kWarning, "Factorization identifies rank deficiency of %d\n", (int)rank_deficiency); return rank_deficiency; } - /** * 2. Elimination other elements by the pivot */ @@ -1311,7 +1295,18 @@ HighsInt HFactor::buildKernel() { HighsInt p2 = p1 + mr_count[iRow]; HighsInt p3 = mr_start[iRow] = mr_index.size(); mr_space[iRow] *= 2; + + HighsInt* mr_index_from_p = mr_index.data(); + mr_index.resize(p3 + mr_space[iRow]); + + HighsInt* mr_index_to_p = mr_index.data(); + if (mr_index_to_p != mr_index_from_p) { + printf("Resize data shift from %p to %p\n", + (void*)mr_index_from_p, (void*)mr_index_to_p); + resize_data_shift = true; + } + copy(&mr_index[p1], &mr_index[p2], &mr_index[p3]); } rowInsert(iCol, iRow); @@ -1343,6 +1338,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; From 31b27b7eb9a41fd077e8c033da728560a2689491 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Sun, 29 Mar 2026 19:01:36 +0100 Subject: [PATCH 14/16] Now to tidy up --- highs/presolve/HPresolve.cpp | 12 ++++++++--- highs/util/HFactor.cpp | 41 +++++++++++++++++++----------------- highs/util/HFactor.h | 2 ++ 3 files changed, 33 insertions(+), 22 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 7f32427744..f86a66daa8 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6519,7 +6519,7 @@ 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 kMaxDependentEquationsTime = 1000; + const double kMaxDependentEquationsTime = 100; const double time_limit = std::max( 1.0, std::min(0.01 * options->time_limit, kMaxDependentEquationsTime)); factor.setTimeLimit(time_limit); @@ -6585,11 +6585,14 @@ HPresolve::Result HPresolve::removeDependentEquations( static_cast(model->num_col_), static_cast(num_nz), static_cast(num_removed_row), static_cast(num_removed_nz), time_taken); + double min_time_bound = kHighsInf; + double max_time_bound = -kHighsInf; + 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 (limit = %.2fs)", + "in %.2fs with bounds in (%.2fs, %.2fs) and limit = %.2fs)", static_cast(num_equations), highsIntToPlural(num_equations).c_str(), static_cast(num_variables), static_cast(model->num_col_), @@ -6597,7 +6600,10 @@ HPresolve::Result HPresolve::removeDependentEquations( num_nz == 1 ? "" : "s", static_cast(num_removed_row), highsIntToPlural(num_removed_row).c_str(), static_cast(num_removed_nz), - highsIntToPlural(num_removed_nz).c_str(), time_taken, time_limit); + highsIntToPlural(num_removed_nz).c_str(), 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", diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index 2e4b67f789..c4fbc65b5d 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -875,6 +875,8 @@ void HFactor::buildSimple() { // 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); } @@ -892,7 +894,7 @@ HighsInt HFactor::buildKernel() { HighsInt num_defer_pivot = 0; // Initial timer frequency: may be reduced if iterations get slow const HighsInt max_timer_frequency = 1000; - HighsInt timer_frequency = 10; + HighsInt timer_frequency = 100; HighsInt log_frequency = 1000; double previous_iteration_time = build_timer_->read(); double average_iteration_time = 0; @@ -901,10 +903,20 @@ HighsInt HFactor::buildKernel() { bool resize_data_shift = false; const bool check_for_timeout = this->time_limit_ < kHighsInf; - auto resizeIntDataShift = [&](std::vector i_vector, const HighsInt to_size) { + auto resizeHighsInt = [&](std::vector& 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; + if (i_vector.data() != from_p) + resize_data_shift = true; + }; + + auto resizeDouble = [&](std::vector& 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; }; const HighsInt check_nwork = -11; @@ -943,16 +955,18 @@ HighsInt HFactor::buildKernel() { HighsInt iterations_left = kernel_dim - search_k + num_defer_pivot + 1; HighsInt iterations_left0 = iterations_left; HighsInt iterations_left1 = 0; - if (average_num_active_nz_change_rate < 0) { + if (average_num_active_nz_change_rate < -1) { HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; iterations_left1 = active_nz_iterations_left; 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; + min_time_bound_ = std::min(total_time_bound, min_time_bound_); + max_time_bound_ = std::max(total_time_bound, max_time_bound_); const bool stop = current_time > this->time_limit_ || total_time_bound > this->time_limit_; - if (search_k % log_frequency == 0 || stop || resize_data_shift) { + if (search_k % log_frequency == 0 || stop || resize_data_shift || total_time_bound < 0) { printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" " ActiveNzRate(%6d, Av = %6d);" " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s%s\n", @@ -1268,8 +1282,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]); @@ -1295,18 +1309,7 @@ HighsInt HFactor::buildKernel() { HighsInt p2 = p1 + mr_count[iRow]; HighsInt p3 = mr_start[iRow] = mr_index.size(); mr_space[iRow] *= 2; - - HighsInt* mr_index_from_p = mr_index.data(); - - mr_index.resize(p3 + mr_space[iRow]); - - HighsInt* mr_index_to_p = mr_index.data(); - if (mr_index_to_p != mr_index_from_p) { - printf("Resize data shift from %p to %p\n", - (void*)mr_index_from_p, (void*)mr_index_to_p); - resize_data_shift = true; - } - + resizeHighsInt(mr_index, p3 + mr_space[iRow]); copy(&mr_index[p1], &mr_index[p2], &mr_index[p3]); } rowInsert(iCol, iRow); diff --git a/highs/util/HFactor.h b/highs/util/HFactor.h index 7cd7e2d78d..1ef4408477 100644 --- a/highs/util/HFactor.h +++ b/highs/util/HFactor.h @@ -337,6 +337,8 @@ class HFactor { HighsInt kernel_dim; HighsInt kernel_num_el; HighsInt num_active_nz_; + double min_time_bound_; + double max_time_bound_; /** * Data of the factor From 1ebb58146284c7ae7b8484ce6261ad68d2551435 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Sun, 29 Mar 2026 22:36:24 +0100 Subject: [PATCH 15/16] Formatted --- highs/presolve/HPresolve.cpp | 28 +++--- highs/simplex/HighsSimplexAnalysis.cpp | 4 +- highs/util/HFactor.cpp | 134 ++++++++++++++----------- highs/util/HFactor.h | 2 +- 4 files changed, 92 insertions(+), 76 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index f86a66daa8..7e0ef5cab1 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6579,15 +6579,15 @@ HPresolve::Result HPresolve::removeDependentEquations( } if (!silent) { highsLogDev(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(num_equations), - static_cast(num_variables), - static_cast(model->num_col_), static_cast(num_nz), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); - double min_time_bound = kHighsInf; - double max_time_bound = -kHighsInf; - + "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", + model->model_name_.c_str(), static_cast(num_equations), + static_cast(num_variables), + static_cast(model->num_col_), static_cast(num_nz), + static_cast(num_removed_row), + static_cast(num_removed_nz), time_taken); + double min_time_bound = kHighsInf; + double max_time_bound = -kHighsInf; + highsLogUser( options->log_options, HighsLogType::kInfo, "Search of %d equation%s with %d / %d variable%s and %d nonzero%s " @@ -6601,14 +6601,12 @@ HPresolve::Result HPresolve::removeDependentEquations( highsIntToPlural(num_removed_row).c_str(), static_cast(num_removed_nz), highsIntToPlural(num_removed_nz).c_str(), time_taken, - factor.min_time_bound_, - factor.max_time_bound_, - time_limit); + 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(num_fictitious_rows_skipped), - highsIntToPlural(num_fictitious_rows_skipped).c_str()); + ", avoiding %d fictitious row%s", + static_cast(num_fictitious_rows_skipped), + highsIntToPlural(num_fictitious_rows_skipped).c_str()); highsLogUser(options->log_options, HighsLogType::kInfo, "\n"); } return returnOk(); diff --git a/highs/simplex/HighsSimplexAnalysis.cpp b/highs/simplex/HighsSimplexAnalysis.cpp index dcaf2ae545..1e9f1aba80 100644 --- a/highs/simplex/HighsSimplexAnalysis.cpp +++ b/highs/simplex/HighsSimplexAnalysis.cpp @@ -1229,8 +1229,8 @@ void HighsSimplexAnalysis::updateInvertFormData(const HFactor& factor) { factor.invert_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 + factor.kernel_dim); + 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 c4fbc65b5d..e5617cb626 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -823,7 +823,8 @@ 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 + // 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++) { @@ -860,11 +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 + // In the active part of the kernel colInsert(iCol, iRow, value); rowInsert(iCol, iRow); } else { - // Above the active part of the kernel + // Above the active part of the kernel colStoreN(iCol, iRow, value); } } @@ -872,7 +873,7 @@ void HFactor::buildSimple() { clinkAdd(iCol, mc_count_a[iCol]); } build_synthetic_tick += (num_row + nwork + MCcountX) * 40 + mr_countX * 20; - // Record the dimension and number of entries in the kernel + // Record the dimension and number of entries in the kernel. kernel_dim = nwork; kernel_num_el = num_active_nz_; min_time_bound_ = kHighsInf; @@ -887,108 +888,119 @@ 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; - HighsInt search_k = 0; + // 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 + 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; - HighsInt log_frequency = 1000; + // Need to maintain an averate iteration time double previous_iteration_time = build_timer_->read(); double average_iteration_time = 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; - const bool check_for_timeout = this->time_limit_ < kHighsInf; + // 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& i_vector, - const HighsInt to_size) { + 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; + if (i_vector.data() != from_p) resize_data_shift = true; }; - + auto resizeDouble = [&](std::vector& d_vector, - const HighsInt to_size) { + 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; + 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 + // 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); - const double mu0 = (1.0 * timer_frequency) / (10.0 * max_timer_frequency); - assert(mu0 <= 0.2); - const double mu1 = 1.0 - mu0; if (!resize_data_shift) { - average_iteration_time = - mu1 * average_iteration_time + mu0 * iteration_time; - - // Determine whether to reduce the timer frequency - if (time_difference > this->time_limit_ / 1e1) - timer_frequency = std::max(HighsInt(1), timer_frequency / 10); + 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(num_active_nz_ - previous_num_active_nz) / - static_cast(timer_frequency); + static_cast(num_active_nz_ - previous_num_active_nz) / + static_cast(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; + 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; - HighsInt iterations_left0 = iterations_left; - HighsInt iterations_left1 = 0; if (average_num_active_nz_change_rate < -1) { - HighsInt active_nz_iterations_left = -num_active_nz_ / average_num_active_nz_change_rate; - iterations_left1 = active_nz_iterations_left; - iterations_left = std::min(active_nz_iterations_left, iterations_left); + 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_); - const bool stop = current_time > this->time_limit_ || - total_time_bound > this->time_limit_; - if (search_k % log_frequency == 0 || stop || resize_data_shift || total_time_bound < 0) { - printf("HFactor::buildKernel k = %7d; Kdim = %7d; Defer =%4d; left(%7d, %7d, %7d); num_active_nz_ = %7d;" - " ActiveNzRate(%6d, Av = %6d);" - " Time (Cu = %6.2f; IterMs(%.6f; Av = %.6f); Bd = %6.2f Lim = %6.2f) Fq = %4d%s%s\n", - int(search_k), int(kernel_dim), - int(num_defer_pivot), - int(iterations_left0), - int(iterations_left1), - int(iterations_left), - int(num_active_nz_), - int(num_active_nz_change_rate), - int(average_num_active_nz_change_rate), - current_time, 1000*iteration_time, 1000*average_iteration_time, - total_time_bound, this->time_limit_, - int(timer_frequency), - resize_data_shift ? " - ReZ-shift " : "", - stop ? " - Stop" : ""); - } - /* + // 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; } /** @@ -1142,6 +1154,8 @@ 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 @@ -1205,8 +1219,9 @@ HighsInt HFactor::buildKernel() { l_value.push_back(value); mr_count_before[iRow] = mr_count[iRow]; rowDelete(jColPivot, (int)iRow); - num_active_nz_--; } + // 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]; @@ -1233,6 +1248,7 @@ 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 @@ -1252,6 +1268,7 @@ 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; @@ -1295,8 +1312,9 @@ HighsInt HFactor::buildKernel() { HighsInt iRow = mwz_column_index[i]; if (mwz_column_mark[iRow]) { colInsert(iCol, iRow, -my_pivot * mwz_column_array[iRow]); - num_active_nz_++; - } + // One active entry is gained for each instance of fill-in + num_active_nz_++; + } } // 2.4.4.3 Fill into the row copy diff --git a/highs/util/HFactor.h b/highs/util/HFactor.h index 1ef4408477..3dfae4cfab 100644 --- a/highs/util/HFactor.h +++ b/highs/util/HFactor.h @@ -339,7 +339,7 @@ class HFactor { HighsInt num_active_nz_; double min_time_bound_; double max_time_bound_; - + /** * Data of the factor */ From 25e2870186c1bc2c887eb6e7de1f4a10ceb0ba73 Mon Sep 17 00:00:00 2001 From: Julian Hall Date: Sun, 29 Mar 2026 23:08:44 +0100 Subject: [PATCH 16/16] Cleaned up dev code --- highs/presolve/HPresolve.cpp | 42 +++++++++--------------------------- highs/util/HFactor.cpp | 3 +++ 2 files changed, 13 insertions(+), 32 deletions(-) diff --git a/highs/presolve/HPresolve.cpp b/highs/presolve/HPresolve.cpp index 7e0ef5cab1..12e5e94e60 100644 --- a/highs/presolve/HPresolve.cpp +++ b/highs/presolve/HPresolve.cpp @@ -6485,26 +6485,14 @@ HPresolve::Result HPresolve::removeDependentEquations( for (HighsInt iCol = 0; iCol < model->num_col_; iCol++) if (row_count[iCol]) num_variables++; HighsInt num_nz = matrix.numNz(); - // - // Although two duplicated equations in any number of variables lead - // to dependency, it's not worth looking for dependent equations if - // there are few relative to the number of variables - // - // The following heuristic picks up all Mittelmann problems with - // meaningful reductions - const bool not_consider_dependent_equations = - num_equations < 0.5 * num_variables; const bool silent = silentLog(); if (!silent) highsLogUser(options->log_options, HighsLogType::kInfo, - "%sonsidering dependency of %d equation%s in %d variable%s " + "Considering dependency of %d equation%s in %d variable%s " "with %d nonzero%s\n", - not_consider_dependent_equations ? "Not c" : "C", int(num_equations), highsIntToPlural(num_equations).c_str(), int(num_variables), highsIntToPlural(num_variables).c_str(), int(num_nz), highsIntToPlural(num_nz).c_str()); - // if (not_consider_dependent_equations) return returnOk(); - // // Identify any dependent equations std::vector colSet(num_equations); std::iota(colSet.begin(), colSet.end(), 0); @@ -6578,30 +6566,20 @@ HPresolve::Result HPresolve::removeDependentEquations( } } if (!silent) { - highsLogDev(options->log_options, HighsLogType::kInfo, - "GrepDependentEq,%s,%d,%d,%d,%d,%d,%d,%g\n", - model->model_name_.c_str(), static_cast(num_equations), - static_cast(num_variables), - static_cast(model->num_col_), static_cast(num_nz), - static_cast(num_removed_row), - static_cast(num_removed_nz), time_taken); - double min_time_bound = kHighsInf; - double max_time_bound = -kHighsInf; - 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)", - static_cast(num_equations), - highsIntToPlural(num_equations).c_str(), - static_cast(num_variables), static_cast(model->num_col_), - highsIntToPlural(num_variables).c_str(), static_cast(num_nz), - num_nz == 1 ? "" : "s", static_cast(num_removed_row), - highsIntToPlural(num_removed_row).c_str(), - static_cast(num_removed_nz), - highsIntToPlural(num_removed_nz).c_str(), time_taken, - factor.min_time_bound_, factor.max_time_bound_, time_limit); + // clang-format off + static_cast(num_equations), highsIntToPlural(num_equations).c_str(), + static_cast(num_variables), + static_cast(model->num_col_), highsIntToPlural(num_variables).c_str(), + static_cast(num_nz), highsIntToPlural(num_nz).c_str(), + static_cast(num_removed_row), highsIntToPlural(num_removed_row).c_str(), + static_cast(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", diff --git a/highs/util/HFactor.cpp b/highs/util/HFactor.cpp index e5617cb626..38b02ccc8d 100644 --- a/highs/util/HFactor.cpp +++ b/highs/util/HFactor.cpp @@ -896,6 +896,9 @@ HighsInt HFactor::buildKernel() { // 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