merge hpm updates

Author: galabovaa

docs/src/options/definitions.md

@@ -6,7 +6,7 @@
 - Default: "choose"
 
 ## [solver](@id option-solver)
-- Solver option: "simplex", "choose", "ipm" or "pdlp". If "simplex"/"ipm"/"pdlp" is chosen then, for a MIP (QP) the integrality constraint (quadratic term) will be ignored
+- Solver option: "simplex", "choose", "ipm", "hipo" or "pdlp". If "simplex"/"ipm"/"hipo"/"pdlp" is chosen then, for a MIP (QP) the integrality constraint (quadratic term) will be ignored
 - Type: string
 - Default: "choose"
 
@@ -409,6 +409,27 @@
 - Range: {0, 2147483647}
 - Default: 2147483647
 
+## [hipo\_system](@id option-hipo-system)
+- Type of Newton system for HiPO: "augmented", "normaleq" or "choose"
+- Type: string
+- Default: "choose"
+
+## [hipo\_parallel\_type](@id option-hipo-parallel)
+- Type of parallelism for HiPO: "tree", "node", "both"
+- Type: string
+- Default: "both"
+
+## hipo\_block\_size
+- Block size for dense linear algebra within HiPO
+- Type: integer
+- Range: {0, 2147483647}
+- Default: 128
+
+## pdlp\_native\_termination
+- Use native termination for PDLP solver: Default = false
+- Type: boolean
+- Default: "false"
+
 ## pdlp\_scaling
 - Scaling option for PDLP solver: Default = true
 - Type: boolean

docs/src/parallel.md

@@ -4,10 +4,10 @@
 
 HiGHS currently has limited opportunities for exploiting parallel
 computing. When using a CPU, these are currently restricted to the
-dual simplex solver for LP, and the MIP solver. Details of these and
-future plans are set out below. HiGHS has an implementation of a first
-order method (PDLP) for solving LPs that can exploit the availability of a
-[GPU](@ref gpu).
+dual simplex solver for LP, the factorisation-based interior point solver,
+and the MIP solver. Details of these and future plans are set out below. 
+HiGHS has an implementation of a first order method (PDLP) for solving LPs 
+that can exploit the availability of a [GPU](@ref highs-gpu).
 
 By default, when running in parallel, HiGHS will use half the
 available threads on a machine. This number can be modified by setting
@@ -47,14 +47,32 @@ clique tables, and when the interior point solver is used to compute
 the analytic centre. This parallelism is always advantageous, so is
 performed regardless of the value of the [parallel](@ref) option.
 
+## IPM
+
+The interior point solver HiPO uses multiple threads to process the 
+elimination tree during the multifrontal factorisation (_tree level_)
+and to perform the dense factorisation of the frontal matrices 
+(_node level_).
+
+If the [parallel](@ref) option is set "on", the level of parallelism is 
+determined by the [hipo\_parallel\_type](@ref option-hipo-parallel) option, 
+which can be "tree" for tree level only, "node" for node level only, or 
+"both" for both levels.
+
+If the [parallel](@ref) option is set "choose", the solver selects which 
+level to use based on a heuristic. When the [parallel](@ref) option is set 
+"choose" or "off", the value of the hipo\_parallel\_type option is ignored.
+
+
 ## Future plans
 
 The MIP solver has been written with parallel tree search in mind. Some
 work has started (Feb 2025), and it is hoped that a prototype solver
 will be available during 2025.
 
-Development of a new parallel interior point solver began 2023, and a
-prototype solver is expected to be be available during 2025.
+Multi-threading within HiPO will be extended to other phases of the solver,
+including the solve phase of the factorisation and the process of assemblying 
+the matrices.
 
 First-order solvers for LP are still very much in their infancy, and
 are not robust. Hence the availability of a PDLP solver for LP is

docs/src/solvers.md

@@ -26,15 +26,30 @@ J. A. J. Hall, Mathematical Programming Computation, 10 (1), 119-142,
 
 #### Interior point
 
-HiGHS has one interior point (IPM) solver based on the preconditioned conjugate gradient method, as discussed in
+HiGHS has two interior point (IPM) solvers:
 
-_Implementation of an interior point method with basis
-preconditioning_, L. Schork and J. Gondzio, Mathematical Programming Computation, 12, 603-635, 2020. [DOI:
-10.1007/s12532-020-00181-8](https://link.springer.com/article/10.1007/s12532-020-00181-8).
+* IPX is based on the preconditioned conjugate gradient method, as discussed in
 
-This solver is serial. An interior point solver based on direct factorization is being developed.
+  _Implementation of an interior point method with basis
+  preconditioning_, Mathematical Programming Computation, 12, 603-635, 2020. [DOI:
+  10.1007/s12532-020-00181-8](https://link.springer.com/article/10.1007/s12532-020-00181-8).
 
-Setting the option [__solver__](@ref option-solver) to "ipm" forces the IPM solver to be used
+  This solver is serial.
+
+  Setting the option [__solver__](@ref option-solver) to "ipm" forces the IPX solver to be used
+
+* HiPO is based on a direct factorisation, as discussed in 
+
+  _A factorisation-based regularised interior point method using the augmented system_, F. Zanetti and J. Gondzio, 2025, 
+  [available on arxiv](https://arxiv.org/abs/2508.04370)
+
+  This solver is parallel.
+
+  Setting the option [__solver__](@ref option-solver) to "hipo" forces the HiPO solver to be used.
+
+  The [hipo\_system](@ref option-hipo-system) option can be used to select the approach to use when solving the Newton systems 
+  within the interior point solver: select "augmented" to force the solver to use the augmented system, "normaleq" for normal 
+  equations, or "choose" to leave the choice to the solver.
 
 #### Primal-dual hybrid gradient method
 
@@ -55,7 +70,7 @@ search is work in progress.
 ## QP
 
 The HiGHS solver for convex QP problems uses an established primal
-active set method. The new interior point solver will also be able to
+active set method. The new interior point solver HiPO will soon be able to
 solve convex QP problems.
 
 

highs/ipm/IpxWrapper.cpp

@@ -499,11 +499,22 @@ HighsStatus solveLpHipo(const HighsOptions& options, HighsTimer& timer,
   // hipo_options.refine_with_ipx = true;
   // hipo_options.display_ipx = true;
 
-  // Option parallel for now is just "on", "off", "choose".
-  // hipo can accept also partially on, to select only tree or node parallelism.
-  // It is worth considering whether this choice should be exposed to the user.
-  if (options.parallel == kHighsOnString)
-    hipo_options.parallel = hipo::kOptionParallelOn;
+  // if option parallel is on, it can be refined by option hipo_parallel_type
+  if (options.parallel == kHighsOnString) {
+    if (options.hipo_parallel_type == kHipoTreeString)
+      hipo_options.parallel = hipo::kOptionParallelTreeOnly;
+    else if (options.hipo_parallel_type == kHipoNodeString)
+      hipo_options.parallel = hipo::kOptionParallelNodeOnly;
+    else if (options.hipo_parallel_type == kHipoBothString)
+      hipo_options.parallel = hipo::kOptionParallelOn;
+    else {
+      highsLogUser(options.log_options, HighsLogType::kError,
+                   "Unknown value of option %s\n", kHipoParallelString.c_str());
+      model_status = HighsModelStatus::kSolveError;
+      return HighsStatus::kError;
+    }
+  }
+  // otherwise, option hipo_parallel_type is ignored
   else if (options.parallel == kHighsOffString)
     hipo_options.parallel = hipo::kOptionParallelOff;
   else {
@@ -520,18 +531,13 @@ HighsStatus solveLpHipo(const HighsOptions& options, HighsTimer& timer,
     hipo_options.nla = hipo::kOptionNlaChoose;
   } else {
     highsLogUser(options.log_options, HighsLogType::kError,
-                 "Unknown value of option hipo_system\n");
+                 "Unknown value of option %s\n", kHipoSystemString.c_str());
     model_status = HighsModelStatus::kSolveError;
     return HighsStatus::kError;
   }
 
-  // ===========================================================================
-  // TO DO
-  // - consider adding options for parallel tree/node
-  // - block size for dense factorisation can have large impact on performance
-  //   and depends on the specific architecture. It may be worth exposing it to
-  //   the user as an advanced option.
-  // ===========================================================================
+  // block size option
+  hipo_options.block_size = options.hipo_block_size;
 
   hipo.set(hipo_options, options.log_options, callback, timer);
 

highs/ipm/hipo/factorhighs/Analyse.cpp

@@ -17,7 +17,7 @@
 namespace hipo {
 
 Analyse::Analyse(const std::vector<Int>& rows, const std::vector<Int>& ptr,
-                 const std::vector<Int>& signs, const Log* log,
+                 const std::vector<Int>& signs, Int nb, const Log* log,
                  DataCollector& data)
     : log_{log}, data_{data} {
   // Input the symmetric matrix to be analysed in CSC format.
@@ -29,7 +29,7 @@ Analyse::Analyse(const std::vector<Int>& rows, const std::vector<Int>& ptr,
   n_ = ptr.size() - 1;
   nz_ = rows.size();
   signs_ = signs;
-  nb_ = kBlockSize;
+  nb_ = nb;
 
   // Create upper triangular part
   rows_upper_.resize(nz_);
@@ -111,12 +111,12 @@ Int Analyse::getPermutation() {
   // fix seed of rng inside Metis, to make it deterministic (?)
   options[METIS_OPTION_SEED] = 42;
 
-  if (log_) log_->printDevInfo("Metis...");
+  if (log_) log_->printDevInfo("Metis...\n");
 
   Int status = METIS_NodeND(&n_, temp_ptr.data(), temp_rows.data(), NULL,
                             options, perm_.data(), iperm_.data());
 
-  if (log_) log_->printDevInfo("done\n");
+  if (log_) log_->printDevInfo("...done\n");
   if (status != METIS_OK) {
     if (log_) log_->printDevInfo("Error with Metis\n");
     return kRetMetisError;
@@ -1357,18 +1357,9 @@ Int Analyse::run(Symbolic& S) {
   computeBlockStart();
   computeCriticalPath();
 
-  // Too many nonzeros for the integer type selected
-  if (nz_factor_ >= std::numeric_limits<Int>::max()) {
-    if (log_) log_->printDevInfo("Integer overflow in analyse phase\n");
-    return kRetIntOverflow;
-  }
-
   // move relevant stuff into S
   S.n_ = n_;
   S.sn_ = sn_count_;
-
-  S.sn_ = sn_count_;
-  S.n_ = n_;
   S.nz_ = nz_factor_;
   S.fillin_ = (double)nz_factor_ / nz_;
   S.artificial_nz_ = artificial_nz_;
@@ -1378,6 +1369,7 @@ Int Analyse::run(Symbolic& S) {
   S.largest_front_ = *std::max_element(sn_indices_.begin(), sn_indices_.end());
   S.serial_storage_ = serial_storage_;
   S.flops_ = dense_ops_;
+  S.block_size_ = nb_;
 
   // compute largest supernode
   std::vector<Int> sn_size(sn_start_.begin() + 1, sn_start_.end());
@@ -1391,14 +1383,18 @@ Int Analyse::run(Symbolic& S) {
     if (i <= 100) S.sn_size_100_++;
   }
 
+  // Too many nonzeros for the integer type selected.
+  // Check after statistics have been moved into S, so that info is accessible
+  // for debug logging.
+  if (nz_factor_ >= std::numeric_limits<Int>::max()) {
+    if (log_) log_->printDevInfo("Integer overflow in analyse phase\n");
+    return kRetIntOverflow;
+  }
+
   // permute signs of pivots
   S.pivot_sign_ = std::move(signs_);
   permuteVector(S.pivot_sign_, perm_);
 
-  S.nz_ = nz_factor_;
-  S.flops_ = dense_ops_;
-  S.spops_ = sparse_ops_;
-  S.critops_ = critical_ops_;
   S.iperm_ = std::move(iperm_);
   S.rows_ = std::move(rows_sn_);
   S.ptr_ = std::move(ptr_sn_);

highs/ipm/hipo/factorhighs/Analyse.h

@@ -106,7 +106,8 @@ class Analyse {
  public:
   // Constructor: matrix must be in lower triangular format
   Analyse(const std::vector<Int>& rows, const std::vector<Int>& ptr,
-          const std::vector<Int>& signs, const Log* log, DataCollector& data);
+          const std::vector<Int>& signs, Int nb, const Log* log,
+          DataCollector& data);
 
   // Run analyse phase and save the result in Symbolic object S
   Int run(Symbolic& S);

highs/ipm/hipo/factorhighs/FactorHiGHS.cpp

@@ -20,8 +20,10 @@ FHsolver::FHsolver(const Log* log) : log_{log} {
 }
 
 FHsolver::~FHsolver() {
-  data_.printTimes(*log_);
-  data_.printIter(*log_);
+  if (log_) {
+    data_.printTimes(*log_);
+    data_.printIter(*log_);
+  }
 }
 
 void FHsolver::newIter() { data_.append(); }
@@ -31,10 +33,12 @@ void FHsolver::setRegularisation(double reg_p, double reg_d) {
   regul_.dual = reg_d;
 }
 
+void FHsolver::setBlockSize(Int nb) { nb_ = nb; }
+
 Int FHsolver::analyse(Symbolic& S, const std::vector<Int>& rows,
                       const std::vector<Int>& ptr,
                       const std::vector<Int>& signs) {
-  Analyse an_obj(rows, ptr, signs, log_, data_);
+  Analyse an_obj(rows, ptr, signs, nb_, log_, data_);
   return an_obj.run(S);
 }
 

highs/ipm/hipo/factorhighs/FactorHiGHS.h

@@ -52,6 +52,10 @@ To add static regularisation when the pivots are selected, use
 setRegularisation(reg_p,reg_d) to choose values of primal and dual
 regularisation. If regularisation is already added to the matrix, ignore.
 
+The default block size is 128. To set a different block size, use setBlockSize.
+Make sure that the block size used for the analyse and factorise phases are the
+same.
+
 */
 
 namespace hipo {
@@ -61,6 +65,8 @@ class FHsolver {
   DataCollector data_;
   Regul regul_;
 
+  Int nb_ = 128;  // block size
+
  public:
   // Create object and initialise DataCollector
   FHsolver(const Log* log = nullptr);
@@ -89,6 +95,9 @@ class FHsolver {
   // Set values for static regularisation to be added when a pivot is selected.
   // If regularisation is already added to the matrix, ignore.
   void setRegularisation(double reg_p, double reg_d);
+
+  // Set the block size for dense linear algebra.
+  void setBlockSize(Int nb);
 };
 
 }  // namespace hipo

highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h

@@ -40,7 +40,6 @@ const Int kMaxIterRelax = 10;
 const Int kSnSizeRelax = 16;
 
 // dense factorisation
-const Int kBlockSize = 128;
 const double kAlphaBK = 0.01;  //(sqrt(17.0) + 1.0) / 8.0;
 const Int kBlockGrainSize = 1;
 const Int kBlockParallelThreshold = 5;

highs/ipm/hipo/factorhighs/Symbolic.cpp

@@ -7,7 +7,7 @@
 
 namespace hipo {
 
-Symbolic::Symbolic() : block_size_{kBlockSize} {}
+Symbolic::Symbolic() {}
 
 void Symbolic::setParallel(bool par_tree, bool par_node) {
   parallel_tree_ = par_tree;
@@ -71,7 +71,7 @@ void Symbolic::print(const Log& log, bool verbose) const {
   log_stream << textline("Flops:") << sci(flops_, 0, 1) << '\n';
   if (verbose) {
     log_stream << textline("Sparse ops:") << sci(spops_, 0, 1) << '\n';
-    log_stream << textline("critical ops:") << sci(critops_, 0, 1) << '\n';
+    log_stream << textline("Critical ops:") << sci(critops_, 0, 1) << '\n';
     log_stream << textline("Max tree speedup:") << fix(flops_ / critops_, 0, 2)
                << '\n';
     log_stream << textline("Artificial nz:") << sci(artificial_nz_, 0, 1)
@@ -87,7 +87,8 @@ void Symbolic::print(const Log& log, bool verbose) const {
     log_stream << textline("Sn size <= 10:") << integer(sn_size_10_, 0) << '\n';
     log_stream << textline("Sn size <= 100:") << integer(sn_size_100_, 0)
                << '\n';
-    log_stream << textline("Sn avg size:") << sci(n_, 0, 1) << '\n';
+    log_stream << textline("Sn avg size:") << sci((double)n_ / sn_, 0, 1)
+               << '\n';
   }
 
   log_stream << '\n';