From 304ece666b91f27a8d17265812886d472fe25100 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 10 Mar 2026 10:22:00 +0000 Subject: [PATCH 01/23] Move permute to auxiliary --- highs/ipm/hipo/auxiliary/Auxiliary.cpp | 71 ++++++++++++++++++++++++ highs/ipm/hipo/auxiliary/Auxiliary.h | 2 + highs/ipm/hipo/factorhighs/Analyse.cpp | 63 +-------------------- highs/ipm/hipo/factorhighs/Factorise.cpp | 65 +--------------------- 4 files changed, 76 insertions(+), 125 deletions(-) diff --git a/highs/ipm/hipo/auxiliary/Auxiliary.cpp b/highs/ipm/hipo/auxiliary/Auxiliary.cpp index 547a498c29..3d2a004e98 100644 --- a/highs/ipm/hipo/auxiliary/Auxiliary.cpp +++ b/highs/ipm/hipo/auxiliary/Auxiliary.cpp @@ -82,6 +82,77 @@ void transpose(const std::vector& ptr, const std::vector& rows, } } +void permuteSym(const std::vector& iperm, std::vector& ptr, + std::vector& rows, std::vector& val, bool lower) { + // Symmetric permutation of the lower (upper) triangular matrix M based on + // inverse permutation iperm. The resulting matrix is lower (upper) triangular, + // regardless of the input matrix. + + const Int n = ptr.size() - 1; + std::vector work(n, 0); + const bool use_val = !val.empty(); + + // go through the columns to count the nonzeros + for (Int j = 0; j < n; ++j) { + // get new index of column + const Int col = iperm[j]; + + // go through elements of column + for (Int el = ptr[j]; el < ptr[j + 1]; ++el) { + const Int i = rows[el]; + + // ignore potential entries in upper(lower) triangular part + if ((lower && i < j) || (!lower && i > j)) continue; + + // get new index of row + const Int row = iperm[i]; + + // if only lower triangular part is used, col is smaller than row + Int actual_col = lower ? std::min(row, col) : std::max(row, col); + ++work[actual_col]; + } + } + + std::vector new_ptr(n + 1); + + // get column pointers by summing the count of nonzeros in each column. + // copy column pointers into work + counts2Ptr(new_ptr, work); + + std::vector new_rows(new_ptr.back()); + std::vector new_val; + if (use_val) new_val.resize(new_ptr.back()); + + // go through the columns to assign row indices + for (Int j = 0; j < n; ++j) { + // get new index of column + const Int col = iperm[j]; + + // go through elements of column + for (Int el = ptr[j]; el < ptr[j + 1]; ++el) { + const Int i = rows[el]; + + // ignore potential entries in upper triangular part + if ((lower && i < j) || (!lower && i > j)) continue; + + // get new index of row + const Int row = iperm[i]; + + // if only lower triangular part is used, col is smaller than row + const Int actual_col = lower ? std::min(row, col) : std::max(row, col); + const Int actual_row = lower ? std::max(row, col) : std::min(row, col); + + Int pos = work[actual_col]++; + new_rows[pos] = actual_row; + if (use_val) new_val[pos] = val[el]; + } + } + + ptr = std::move(new_ptr); + rows = std::move(new_rows); + if (use_val) val = std::move(new_val); +} + void childrenLinkedList(const std::vector& parent, std::vector& head, std::vector& next) { // Create linked lists of children in elimination tree. diff --git a/highs/ipm/hipo/auxiliary/Auxiliary.h b/highs/ipm/hipo/auxiliary/Auxiliary.h index b5ea96b6bf..2a7b94f5f7 100644 --- a/highs/ipm/hipo/auxiliary/Auxiliary.h +++ b/highs/ipm/hipo/auxiliary/Auxiliary.h @@ -19,6 +19,8 @@ void transpose(const std::vector& ptr, const std::vector& rows, void transpose(const std::vector& ptr, const std::vector& rows, const std::vector& val, std::vector& ptrT, std::vector& rowsT, std::vector& valT); +void permuteSym(const std::vector& iperm, std::vector& ptr, + std::vector& rows, std::vector& val, bool lower); void childrenLinkedList(const std::vector& parent, std::vector& head, std::vector& next); void reverseLinkedList(std::vector& head, std::vector& next); diff --git a/highs/ipm/hipo/factorhighs/Analyse.cpp b/highs/ipm/hipo/factorhighs/Analyse.cpp index 587c71a845..f4153e72a2 100644 --- a/highs/ipm/hipo/factorhighs/Analyse.cpp +++ b/highs/ipm/hipo/factorhighs/Analyse.cpp @@ -67,67 +67,8 @@ Analyse::Analyse(const std::vector& rows, const std::vector& ptr, } void Analyse::permute(const std::vector& iperm) { - // Symmetric permutation of the upper triangular matrix based on inverse - // permutation iperm. - // The resulting matrix is upper triangular, regardless of the input matrix. - - std::vector work(n_, 0); - - // go through the columns to count the nonzeros - for (Int j = 0; j < n_; ++j) { - // get new index of column - const Int col = iperm[j]; - - // go through elements of column - for (Int el = ptr_upper_[j]; el < ptr_upper_[j + 1]; ++el) { - const Int i = rows_upper_[el]; - - // ignore potential entries in lower triangular part - if (i > j) continue; - - // get new index of row - const Int row = iperm[i]; - - // since only upper triangular part is used, col is larger than row - Int actual_col = std::max(row, col); - ++work[actual_col]; - } - } - - std::vector new_ptr(n_ + 1); - - // get column pointers by summing the count of nonzeros in each column. - // copy column pointers into work - counts2Ptr(new_ptr, work); - - std::vector new_rows(new_ptr.back()); - - // go through the columns to assign row indices - for (Int j = 0; j < n_; ++j) { - // get new index of column - const Int col = iperm[j]; - - // go through elements of column - for (Int el = ptr_upper_[j]; el < ptr_upper_[j + 1]; ++el) { - const Int i = rows_upper_[el]; - - // ignore potential entries in lower triangular part - if (i > j) continue; - - // get new index of row - const Int row = iperm[i]; - - // since only upper triangular part is used, column is larger than row - const Int actual_col = std::max(row, col); - const Int actual_row = std::min(row, col); - - Int pos = work[actual_col]++; - new_rows[pos] = actual_row; - } - } - - ptr_upper_ = std::move(new_ptr); - rows_upper_ = std::move(new_rows); + std::vector empty_vals; + permuteSym(iperm, ptr_upper_, rows_upper_, empty_vals, false); } void Analyse::eTree() { diff --git a/highs/ipm/hipo/factorhighs/Factorise.cpp b/highs/ipm/hipo/factorhighs/Factorise.cpp index af06414b82..d03e203d84 100644 --- a/highs/ipm/hipo/factorhighs/Factorise.cpp +++ b/highs/ipm/hipo/factorhighs/Factorise.cpp @@ -91,70 +91,7 @@ Factorise::Factorise(const Symbolic& S, const std::vector& rowsM, } void Factorise::permute(const std::vector& iperm) { - // Symmetric permutation of the lower triangular matrix M based on inverse - // permutation iperm. - // The resulting matrix is lower triangular, regardless of the input matrix. - - std::vector work(n_, 0); - - // go through the columns to count the nonzeros - for (Int j = 0; j < n_; ++j) { - // get new index of column - const Int col = iperm[j]; - - // go through elements of column - for (Int el = ptrM_[j]; el < ptrM_[j + 1]; ++el) { - const Int i = rowsM_[el]; - - // ignore potential entries in upper triangular part - if (i < j) continue; - - // get new index of row - const Int row = iperm[i]; - - // since only lower triangular part is used, col is smaller than row - Int actual_col = std::min(row, col); - ++work[actual_col]; - } - } - - std::vector new_ptr(n_ + 1); - - // get column pointers by summing the count of nonzeros in each column. - // copy column pointers into work - counts2Ptr(new_ptr, work); - - std::vector new_rows(new_ptr.back()); - std::vector new_val(new_ptr.back()); - - // go through the columns to assign row indices - for (Int j = 0; j < n_; ++j) { - // get new index of column - const Int col = iperm[j]; - - // go through elements of column - for (Int el = ptrM_[j]; el < ptrM_[j + 1]; ++el) { - const Int i = rowsM_[el]; - - // ignore potential entries in upper triangular part - if (i < j) continue; - - // get new index of row - const Int row = iperm[i]; - - // since only lower triangular part is used, col is smaller than row - const Int actual_col = std::min(row, col); - const Int actual_row = std::max(row, col); - - Int pos = work[actual_col]++; - new_rows[pos] = actual_row; - new_val[pos] = valM_[el]; - } - } - - ptrM_ = std::move(new_ptr); - rowsM_ = std::move(new_rows); - valM_ = std::move(new_val); + permuteSym(iperm, ptrM_, rowsM_, valM_, true); } void Factorise::processSupernode(Int sn) { From f06000538dd013a0296a9594eaf9cdf65bf34c8b Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 10 Mar 2026 11:04:28 +0000 Subject: [PATCH 02/23] Prepare data structure to extract matrix --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 87 +++++++++++++----------- highs/ipm/hipo/ipm/FactorHiGHSSolver.h | 12 +--- highs/ipm/hipo/ipm/LinearSolver.h | 12 ++++ highs/ipm/hipo/ipm/Solver.cpp | 4 +- highs/ipm/hipo/ipm/Solver.h | 2 + 5 files changed, 66 insertions(+), 51 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index b0ba1d7682..cc5510a61f 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -12,11 +12,13 @@ namespace hipo { -FactorHiGHSSolver::FactorHiGHSSolver(Options& options, const Model& model, +FactorHiGHSSolver::FactorHiGHSSolver(KktMatrix& kkt, Options& options, + const Model& model, const Regularisation& regul, Info& info, IpmData& record, const LogHighs& log) : FH_(&log, options.block_size), S_{}, + kkt_{kkt}, regul_{regul}, info_{info}, data_{record}, @@ -51,42 +53,42 @@ Int FactorHiGHSSolver::buildASstructure() { // AS matrix must fit into HighsInt if ((Int64)nzBlock11 + mA_ + nzA_ >= kHighsIInf) return kStatusOverflow; - ptrAS_.resize(nA_ + mA_ + 1); - rowsAS_.resize(nzBlock11 + nzA_ + mA_); - valAS_.resize(nzBlock11 + nzA_ + mA_); + kkt_.ptrAS.resize(nA_ + mA_ + 1); + kkt_.rowsAS.resize(nzBlock11 + nzA_ + mA_); + kkt_.valAS.resize(nzBlock11 + nzA_ + mA_); Int next = 0; for (Int i = 0; i < nA_; ++i) { // diagonal element - rowsAS_[next] = i; + kkt_.rowsAS[next] = i; next++; // column of Q if (model_.qp()) { assert(Q_.index_[Q_.start_[i]] == i); for (Int el = Q_.start_[i] + 1; el < Q_.start_[i + 1]; ++el) { - rowsAS_[next] = Q_.index_[el]; - valAS_[next] = -Q_.value_[el]; // values of AS that will not change + kkt_.rowsAS[next] = Q_.index_[el]; + kkt_.valAS[next] = -Q_.value_[el]; // values of AS that will not change ++next; } } // column of A for (Int el = A_.start_[i]; el < A_.start_[i + 1]; ++el) { - rowsAS_[next] = nA_ + A_.index_[el]; - valAS_[next] = A_.value_[el]; // values of AS that will not change + kkt_.rowsAS[next] = nA_ + A_.index_[el]; + kkt_.valAS[next] = A_.value_[el]; // values of AS that will not change ++next; } - ptrAS_[i + 1] = next; + kkt_.ptrAS[i + 1] = next; } // 2,2 block for (Int i = 0; i < mA_; ++i) { - rowsAS_[next] = nA_ + i; + kkt_.rowsAS[next] = nA_ + i; ++next; - ptrAS_[nA_ + i + 1] = ptrAS_[nA_ + i] + 1; + kkt_.ptrAS[nA_ + i + 1] = kkt_.ptrAS[nA_ + i] + 1; } info_.AS_structure_time = clock.stop(); @@ -97,11 +99,12 @@ Int FactorHiGHSSolver::buildASstructure() { Int FactorHiGHSSolver::buildASvalues(const std::vector& scaling) { // build AS values that change during iterations. - assert(!ptrAS_.empty() && !rowsAS_.empty()); + assert(!kkt_.ptrAS.empty() && !kkt_.rowsAS.empty()); for (Int i = 0; i < nA_; ++i) { - valAS_[ptrAS_[i]] = scaling.empty() ? -1.0 : -scaling[i]; - if (model_.qp()) valAS_[ptrAS_[i]] -= model_.sense() * model_.Q().diag(i); + kkt_.valAS[kkt_.ptrAS[i]] = scaling.empty() ? -1.0 : -scaling[i]; + if (model_.qp()) + kkt_.valAS[kkt_.ptrAS[i]] -= model_.sense() * model_.Q().diag(i); } return kStatusOk; @@ -142,11 +145,11 @@ Int FactorHiGHSSolver::buildNEstructure() { } } - ptrNE_.clear(); - rowsNE_.clear(); + kkt_.ptrNE.clear(); + kkt_.rowsNE.clear(); // ptr is allocated its exact size - ptrNE_.resize(mA_ + 1, 0); + kkt_.ptrNE.resize(mA_ + 1, 0); // keep track if given entry is nonzero, in column considered std::vector is_nz(mA_, false); @@ -183,18 +186,18 @@ Int FactorHiGHSSolver::buildNEstructure() { // intersection of row with rows below finished. // if the total number of nonzeros exceeds the maximum, return error. - if ((Int64)ptrNE_[row] + nz_in_col >= + if ((Int64)kkt_.ptrNE[row] + nz_in_col >= NE_nz_limit_.load(std::memory_order_relaxed)) return kStatusOverflow; // update pointers - ptrNE_[row + 1] = ptrNE_[row] + nz_in_col; + kkt_.ptrNE[row + 1] = kkt_.ptrNE[row] + nz_in_col; // now assign indices for (Int i = 0; i < nz_in_col; ++i) { Int index = temp_index[i]; // push_back is better then reserve, because the final length is not known - rowsNE_.push_back(index); + kkt_.rowsNE.push_back(index); is_nz[index] = false; } } @@ -207,9 +210,9 @@ Int FactorHiGHSSolver::buildNEstructure() { Int FactorHiGHSSolver::buildNEvalues(const std::vector& scaling) { // given the NE structure already computed, fill in the NE values - assert(!ptrNE_.empty() && !rowsNE_.empty()); + assert(!kkt_.ptrNE.empty() && !kkt_.rowsNE.empty()); - valNE_.resize(rowsNE_.size()); + kkt_.valNE.resize(kkt_.rowsNE.size()); std::vector work(mA_, 0.0); @@ -244,9 +247,9 @@ Int FactorHiGHSSolver::buildNEvalues(const std::vector& scaling) { // intersection of row with rows below finished. // read from work, using indices of column "row" of AAt - for (Int el = ptrNE_[row]; el < ptrNE_[row + 1]; ++el) { - Int index = rowsNE_[el]; - valNE_[el] = work[index]; + for (Int el = kkt_.ptrNE[row]; el < kkt_.ptrNE[row + 1]; ++el) { + Int index = kkt_.rowsNE[el]; + kkt_.valNE[el] = work[index]; work[index] = 0.0; } } @@ -262,7 +265,7 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) { // Perform analyse phase of augmented system and return symbolic factorisation // in object S and the status. - if (rowsAS_.empty() || ptrAS_.empty()) return kStatusErrorAnalyse; + if (kkt_.rowsAS.empty() || kkt_.ptrAS.empty()) return kStatusErrorAnalyse; // create vector of signs of pivots std::vector pivot_signs(nA_ + mA_, -1); @@ -271,8 +274,8 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) { log_.printDevInfo("Performing AS analyse phase\n"); Clock clock; - Int status = - chooseOrdering(rowsAS_, ptrAS_, pivot_signs, S, ordering_AS_, "AS"); + Int status = chooseOrdering(kkt_.rowsAS, kkt_.ptrAS, pivot_signs, S, + ordering_AS_, "AS"); info_.analyse_AS_time += clock.stop(); return status; @@ -283,7 +286,7 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S) { // in object S and the status. Structure of the matrix must be already // computed. - if (rowsNE_.empty() || ptrNE_.empty()) return kStatusErrorAnalyse; + if (kkt_.rowsNE.empty() || kkt_.ptrNE.empty()) return kStatusErrorAnalyse; // create vector of signs of pivots std::vector pivot_signs(mA_, 1); @@ -291,8 +294,8 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S) { log_.printDevInfo("Performing NE analyse phase\n"); Clock clock; - Int status = - chooseOrdering(rowsNE_, ptrNE_, pivot_signs, S, ordering_NE_, "NE"); + Int status = chooseOrdering(kkt_.rowsNE, kkt_.ptrNE, pivot_signs, S, + ordering_NE_, "NE"); info_.analyse_NE_time += clock.stop(); return status; @@ -317,7 +320,8 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { // factorise matrix clock.start(); - if (FH_.factorise(S_, rowsAS_, ptrAS_, valAS_)) return kStatusErrorFactorise; + if (FH_.factorise(S_, kkt_.rowsAS, kkt_.ptrAS, kkt_.valAS)) + return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -340,7 +344,8 @@ Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { // factorise clock.start(); - if (FH_.factorise(S_, rowsNE_, ptrNE_, valNE_)) return kStatusErrorFactorise; + if (FH_.factorise(S_, kkt_.rowsNE, kkt_.ptrNE, kkt_.valNE)) + return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -772,6 +777,8 @@ Int FactorHiGHSSolver::setNla() { << '\n'; log_.print(log_stream); + kkt_.iperm = S_.iperm(); + return kStatusOk; } @@ -869,16 +876,16 @@ void FactorHiGHSSolver::getReg(std::vector& reg) { void FactorHiGHSSolver::freeASmemory() { // Give up memory used for AS. - freeVector(ptrAS_); - freeVector(rowsAS_); - freeVector(valAS_); + freeVector(kkt_.ptrAS); + freeVector(kkt_.rowsAS); + freeVector(kkt_.valAS); } void FactorHiGHSSolver::freeNEmemory() { // Give up memory used for NE. - freeVector(ptrNE_); - freeVector(rowsNE_); - freeVector(valNE_); + freeVector(kkt_.ptrNE); + freeVector(kkt_.rowsNE); + freeVector(kkt_.valNE); freeVector(ptrA_rw_); freeVector(idxA_rw_); freeVector(corr_A_); diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h index 4e2d88328b..c0c6b98561 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h @@ -20,19 +20,13 @@ class FactorHiGHSSolver : public LinearSolver { // symbolic factorisation Symbolic S_; + KktMatrix& kkt_; + // normal equations data - std::vector ptrNE_; - std::vector rowsNE_; - std::vector valNE_; std::vector ptrA_rw_, idxA_rw_; std::vector corr_A_; std::atomic NE_nz_limit_{kHighsIInf}; - // augmented system data - std::vector ptrAS_; - std::vector rowsAS_; - std::vector valAS_; - const Regularisation& regul_; Info& info_; @@ -67,7 +61,7 @@ class FactorHiGHSSolver : public LinearSolver { Int analyseNE(Symbolic& S); public: - FactorHiGHSSolver(Options& options, const Model& model, + FactorHiGHSSolver(KktMatrix& kkt, Options& options, const Model& model, const Regularisation& regul, Info& info, IpmData& record, const LogHighs& log); diff --git a/highs/ipm/hipo/ipm/LinearSolver.h b/highs/ipm/hipo/ipm/LinearSolver.h index f96d5fec6b..a44373fff4 100644 --- a/highs/ipm/hipo/ipm/LinearSolver.h +++ b/highs/ipm/hipo/ipm/LinearSolver.h @@ -76,6 +76,18 @@ class LinearSolver { virtual void getReg(std::vector& reg){}; }; +struct KktMatrix { + std::vector ptrNE; + std::vector rowsNE; + std::vector valNE; + + std::vector ptrAS; + std::vector rowsAS; + std::vector valAS; + + std::vector iperm; +}; + } // namespace hipo #endif \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index f159e042e1..67f5b72c10 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -133,8 +133,8 @@ bool Solver::initialise() { start_time_ = control_.elapsed(); // initialise linear solver - LS_.reset( - new FactorHiGHSSolver(options_, model_, regul_, info_, it_->data, logH_)); + LS_.reset(new FactorHiGHSSolver(kkt_, options_, model_, regul_, info_, + it_->data, logH_)); if (Int status = LS_->setup()) { info_.status = (Status)status; return true; diff --git a/highs/ipm/hipo/ipm/Solver.h b/highs/ipm/hipo/ipm/Solver.h index 7254a957e2..0fcbd7ec13 100644 --- a/highs/ipm/hipo/ipm/Solver.h +++ b/highs/ipm/hipo/ipm/Solver.h @@ -32,6 +32,8 @@ class Solver { // Linear solver interface std::unique_ptr LS_; + KktMatrix kkt_; + // Iterate object interface std::unique_ptr it_; From c8ff047c5a45d4baf90b12fbc3d3933cb63a3dfe Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 10 Mar 2026 12:37:00 +0000 Subject: [PATCH 03/23] Move kkt structure and values out --- cmake/sources.cmake | 2 + highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 274 ++--------------------- highs/ipm/hipo/ipm/FactorHiGHSSolver.h | 22 +- highs/ipm/hipo/ipm/KktMatrix.cpp | 254 +++++++++++++++++++++ highs/ipm/hipo/ipm/KktMatrix.h | 46 ++++ highs/ipm/hipo/ipm/LinearSolver.h | 12 - highs/ipm/hipo/ipm/Solver.cpp | 12 +- highs/ipm/hipo/ipm/Solver.h | 2 +- 8 files changed, 334 insertions(+), 290 deletions(-) create mode 100644 highs/ipm/hipo/ipm/KktMatrix.cpp create mode 100644 highs/ipm/hipo/ipm/KktMatrix.h diff --git a/cmake/sources.cmake b/cmake/sources.cmake index 1e605de4b7..ac872418eb 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -182,6 +182,7 @@ set(hipo_sources ipm/hipo/ipm/FactorHiGHSSolver.cpp ipm/hipo/ipm/Control.cpp ipm/hipo/ipm/Iterate.cpp + ipm/hipo/ipm/KktMatrix.cpp ipm/hipo/ipm/LogHighs.cpp ipm/hipo/ipm/Model.cpp ipm/hipo/ipm/PreProcess.cpp @@ -195,6 +196,7 @@ set(hipo_headers ipm/hipo/ipm/Control.h ipm/hipo/ipm/Info.h ipm/hipo/ipm/Iterate.h + ipm/hipo/ipm/KktMatrix.h ipm/hipo/ipm/LinearSolver.h ipm/hipo/ipm/LogHighs.h ipm/hipo/ipm/Model.h diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index cc5510a61f..4562a2baf3 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -24,12 +24,6 @@ FactorHiGHSSolver::FactorHiGHSSolver(KktMatrix& kkt, Options& options, data_{record}, log_{log}, model_{model}, - A_{model.A()}, - Q_{model.Q()}, - mA_{A_.num_row_}, - nA_{A_.num_col_}, - nzA_{A_.numNz()}, - nzQ_{Q_.numNz()}, options_{options} {} void FactorHiGHSSolver::clear() { @@ -37,226 +31,6 @@ void FactorHiGHSSolver::clear() { FH_.newIter(); } -// ========================================================================= -// Build structure and values of matrices -// ========================================================================= - -Int FactorHiGHSSolver::buildASstructure() { - // Build lower triangular structure of the augmented system. - // Build values of AS that will not change during the iterations. - - Clock clock; - log_.printDevInfo("Building AS structure\n"); - - const Int nzBlock11 = model_.qp() ? nzQ_ : nA_; - - // AS matrix must fit into HighsInt - if ((Int64)nzBlock11 + mA_ + nzA_ >= kHighsIInf) return kStatusOverflow; - - kkt_.ptrAS.resize(nA_ + mA_ + 1); - kkt_.rowsAS.resize(nzBlock11 + nzA_ + mA_); - kkt_.valAS.resize(nzBlock11 + nzA_ + mA_); - - Int next = 0; - - for (Int i = 0; i < nA_; ++i) { - // diagonal element - kkt_.rowsAS[next] = i; - next++; - - // column of Q - if (model_.qp()) { - assert(Q_.index_[Q_.start_[i]] == i); - for (Int el = Q_.start_[i] + 1; el < Q_.start_[i + 1]; ++el) { - kkt_.rowsAS[next] = Q_.index_[el]; - kkt_.valAS[next] = -Q_.value_[el]; // values of AS that will not change - ++next; - } - } - - // column of A - for (Int el = A_.start_[i]; el < A_.start_[i + 1]; ++el) { - kkt_.rowsAS[next] = nA_ + A_.index_[el]; - kkt_.valAS[next] = A_.value_[el]; // values of AS that will not change - ++next; - } - - kkt_.ptrAS[i + 1] = next; - } - - // 2,2 block - for (Int i = 0; i < mA_; ++i) { - kkt_.rowsAS[next] = nA_ + i; - ++next; - kkt_.ptrAS[nA_ + i + 1] = kkt_.ptrAS[nA_ + i] + 1; - } - - info_.AS_structure_time = clock.stop(); - - return kStatusOk; -} - -Int FactorHiGHSSolver::buildASvalues(const std::vector& scaling) { - // build AS values that change during iterations. - - assert(!kkt_.ptrAS.empty() && !kkt_.rowsAS.empty()); - - for (Int i = 0; i < nA_; ++i) { - kkt_.valAS[kkt_.ptrAS[i]] = scaling.empty() ? -1.0 : -scaling[i]; - if (model_.qp()) - kkt_.valAS[kkt_.ptrAS[i]] -= model_.sense() * model_.Q().diag(i); - } - - return kStatusOk; -} - -Int FactorHiGHSSolver::buildNEstructure() { - // Build lower triangular structure of AAt. - // This approach uses a column-wise copy of A, a partial row-wise copy and a - // vector of corresponding indices. - - // NB: A must have sorted columns for this to work - - Clock clock; - log_.printDevInfo("Building NE structure\n"); - - // create partial row-wise representation without values, and array or - // corresponding indices between cw and rw representation - { - ptrA_rw_.assign(mA_ + 1, 0); - idxA_rw_.assign(nzA_, 0); - - // pointers of row-start - for (Int el = 0; el < nzA_; ++el) ptrA_rw_[A_.index_[el] + 1]++; - for (Int i = 0; i < mA_; ++i) ptrA_rw_[i + 1] += ptrA_rw_[i]; - - std::vector temp = ptrA_rw_; - corr_A_.assign(nzA_, 0); - - // rw-indices and corresponding indices created together - for (Int col = 0; col < nA_; ++col) { - for (Int el = A_.start_[col]; el < A_.start_[col + 1]; ++el) { - Int row = A_.index_[el]; - - corr_A_[temp[row]] = el; - idxA_rw_[temp[row]] = col; - temp[row]++; - } - } - } - - kkt_.ptrNE.clear(); - kkt_.rowsNE.clear(); - - // ptr is allocated its exact size - kkt_.ptrNE.resize(mA_ + 1, 0); - - // keep track if given entry is nonzero, in column considered - std::vector is_nz(mA_, false); - - // temporary storage of indices - std::vector temp_index(mA_); - - for (Int row = 0; row < mA_; ++row) { - // go along the entries of the row, and then down each column. - // this builds the lower triangular part of the row-th column of AAt. - - Int nz_in_col = 0; - - for (Int el = ptrA_rw_[row]; el < ptrA_rw_[row + 1]; ++el) { - Int col = idxA_rw_[el]; - Int corr = corr_A_[el]; - - // for each nonzero in the row, go down corresponding column, starting - // from current position - for (Int colEl = corr; colEl < A_.start_[col + 1]; ++colEl) { - Int row2 = A_.index_[colEl]; - - // row2 is guaranteed to be larger or equal than row - // (provided that the columns of A are sorted) - - // save information that there is nonzero in position (row2,row). - if (!is_nz[row2]) { - is_nz[row2] = true; - temp_index[nz_in_col] = row2; - ++nz_in_col; - } - } - } - // intersection of row with rows below finished. - - // if the total number of nonzeros exceeds the maximum, return error. - if ((Int64)kkt_.ptrNE[row] + nz_in_col >= - NE_nz_limit_.load(std::memory_order_relaxed)) - return kStatusOverflow; - - // update pointers - kkt_.ptrNE[row + 1] = kkt_.ptrNE[row] + nz_in_col; - - // now assign indices - for (Int i = 0; i < nz_in_col; ++i) { - Int index = temp_index[i]; - // push_back is better then reserve, because the final length is not known - kkt_.rowsNE.push_back(index); - is_nz[index] = false; - } - } - - info_.NE_structure_time = clock.stop(); - - return kStatusOk; -} - -Int FactorHiGHSSolver::buildNEvalues(const std::vector& scaling) { - // given the NE structure already computed, fill in the NE values - - assert(!kkt_.ptrNE.empty() && !kkt_.rowsNE.empty()); - - kkt_.valNE.resize(kkt_.rowsNE.size()); - - std::vector work(mA_, 0.0); - - for (Int row = 0; row < mA_; ++row) { - // go along the entries of the row, and then down each column. - // this builds the lower triangular part of the row-th column of AAt. - - for (Int el = ptrA_rw_[row]; el < ptrA_rw_[row + 1]; ++el) { - Int col = idxA_rw_[el]; - Int corr = corr_A_[el]; - - double denom = scaling.empty() ? 1.0 : scaling[col]; - denom += regul_.primal; - if (model_.qp()) denom += model_.sense() * model_.Q().diag(col); - - const double mult = 1.0 / denom; - const double row_value = mult * A_.value_[corr]; - - // for each nonzero in the row, go down corresponding column, starting - // from current position - for (Int colEl = corr; colEl < A_.start_[col + 1]; ++colEl) { - Int row2 = A_.index_[colEl]; - - // row2 is guaranteed to be larger or equal than row - // (provided that the columns of A are sorted) - - // compute and accumulate value - double value = row_value * A_.value_[colEl]; - work[row2] += value; - } - } - // intersection of row with rows below finished. - - // read from work, using indices of column "row" of AAt - for (Int el = kkt_.ptrNE[row]; el < kkt_.ptrNE[row + 1]; ++el) { - Int index = kkt_.rowsNE[el]; - kkt_.valNE[el] = work[index]; - work[index] = 0.0; - } - } - - return kStatusOk; -} - // ========================================================================= // Analyse phase // ========================================================================= @@ -267,9 +41,12 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) { if (kkt_.rowsAS.empty() || kkt_.ptrAS.empty()) return kStatusErrorAnalyse; + const Int m = model_.A().num_row_; + const Int n = model_.A().num_col_; + // create vector of signs of pivots - std::vector pivot_signs(nA_ + mA_, -1); - for (Int i = 0; i < mA_; ++i) pivot_signs[nA_ + i] = 1; + std::vector pivot_signs(n + m, -1); + for (Int i = 0; i < m; ++i) pivot_signs[n + i] = 1; log_.printDevInfo("Performing AS analyse phase\n"); @@ -289,7 +66,7 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S) { if (kkt_.rowsNE.empty() || kkt_.ptrNE.empty()) return kStatusErrorAnalyse; // create vector of signs of pivots - std::vector pivot_signs(mA_, 1); + std::vector pivot_signs(model_.A().num_row_, 1); log_.printDevInfo("Performing NE analyse phase\n"); @@ -312,7 +89,7 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { Clock clock; // build matrix - buildASvalues(scaling); + kkt_.buildASvalues(scaling); info_.matrix_time += clock.stop(); // set static regularisation, since it may have changed @@ -336,7 +113,7 @@ Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { Clock clock; // build matrix - buildNEvalues(scaling); + kkt_.buildNEvalues(scaling); info_.matrix_time += clock.stop(); // set static regularisation, since it may have changed @@ -442,7 +219,9 @@ Int FactorHiGHSSolver::chooseNla() { model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; } else { - Int status = buildNEstructure(); + Clock clock; + Int status = kkt_.buildNEstructure(); + info_.NE_structure_time = clock.stop(); if (status) { failure_NE = true; if (status == kStatusOverflow) { @@ -461,7 +240,7 @@ Int FactorHiGHSSolver::chooseNla() { }; auto run_analyse_AS = [&]() { - Int AS_status = buildASstructure(); + Int AS_status = kkt_.buildASstructure(); if (!AS_status) AS_status = analyseAS(symb_AS); if (AS_status) failure_AS = true; if (AS_status == kStatusOverflow) { @@ -473,7 +252,7 @@ Int FactorHiGHSSolver::chooseNla() { // will be preferred, so stop computation of NE. Int64 NE_nz_limit = symb_AS.nz() * kSymbNzMult; if (failure_AS || NE_nz_limit > kHighsIInf) NE_nz_limit = kHighsIInf; - NE_nz_limit_.store(NE_nz_limit, std::memory_order_relaxed); + kkt_.NE_nz_limit.store(NE_nz_limit, std::memory_order_relaxed); }; // In parallel, run AS analyse and build NE structure. NE analyse runs only @@ -540,10 +319,10 @@ Int FactorHiGHSSolver::chooseNla() { if (status == kStatusOk) { if (options_.nla == kHipoAugmentedString) { S_ = std::move(symb_AS); - freeNEmemory(); + kkt_.freeNEmemory(); } else { S_ = std::move(symb_NE); - freeASmemory(); + kkt_.freeASmemory(); } } @@ -575,7 +354,7 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, std::vector failure(orderings_to_try.size(), 0); if (nla == "NE") { - if (ptr.back() >= NE_nz_limit_.load(std::memory_order_relaxed)) { + if (ptr.back() >= kkt_.NE_nz_limit.load(std::memory_order_relaxed)) { log_.printDevInfo("NE interrupted\n"); return kStatusErrorAnalyse; } @@ -741,7 +520,7 @@ Int FactorHiGHSSolver::setNla() { } if (options_.nla == kHipoAugmentedString) { - Int status = buildASstructure(); + Int status = kkt_.buildASstructure(); if (!status) status = analyseAS(S_); if (status == kStatusOverflow) { log_.printe("AS requested, integer overflow\n"); @@ -753,7 +532,7 @@ Int FactorHiGHSSolver::setNla() { log_stream << textline("Newton system:") << "AS requested\n"; } else if (options_.nla == kHipoNormalEqString) { - Int status = buildNEstructure(); + Int status = kkt_.buildNEstructure(); if (!status) status = analyseNE(S_); if (status == kStatusOverflow) { log_.printe("NE requested, integer overflow\n"); @@ -874,21 +653,4 @@ void FactorHiGHSSolver::getReg(std::vector& reg) { FH_.getRegularisation(reg); } -void FactorHiGHSSolver::freeASmemory() { - // Give up memory used for AS. - freeVector(kkt_.ptrAS); - freeVector(kkt_.rowsAS); - freeVector(kkt_.valAS); -} - -void FactorHiGHSSolver::freeNEmemory() { - // Give up memory used for NE. - freeVector(kkt_.ptrNE); - freeVector(kkt_.rowsNE); - freeVector(kkt_.valNE); - freeVector(ptrA_rw_); - freeVector(idxA_rw_); - freeVector(corr_A_); -} - } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h index c0c6b98561..779758b363 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h @@ -6,6 +6,7 @@ #include "Info.h" #include "IpmData.h" +#include "KktMatrix.h" #include "LinearSolver.h" #include "Model.h" #include "ipm/hipo/auxiliary/IntConfig.h" @@ -22,23 +23,13 @@ class FactorHiGHSSolver : public LinearSolver { KktMatrix& kkt_; - // normal equations data - std::vector ptrA_rw_, idxA_rw_; - std::vector corr_A_; - std::atomic NE_nz_limit_{kHighsIInf}; - const Regularisation& regul_; - Info& info_; IpmData& data_; const LogHighs& log_; - const Model& model_; - const HighsSparseMatrix& A_; - const HighsHessian& Q_; - const Int mA_, nA_, nzA_, nzQ_; - Options& options_; + std::string ordering_AS_ = "none"; std::string ordering_NE_ = "none"; @@ -48,15 +39,6 @@ class FactorHiGHSSolver : public LinearSolver { Int chooseOrdering(const std::vector& rows, const std::vector& ptr, const std::vector& signs, Symbolic& S, std::string& ordering, const std::string& nla); - - Int buildNEstructure(); - Int buildNEvalues(const std::vector& scaling); - void freeNEmemory(); - - Int buildASstructure(); - Int buildASvalues(const std::vector& scaling); - void freeASmemory(); - Int analyseAS(Symbolic& S); Int analyseNE(Symbolic& S); diff --git a/highs/ipm/hipo/ipm/KktMatrix.cpp b/highs/ipm/hipo/ipm/KktMatrix.cpp new file mode 100644 index 0000000000..8dc050ed56 --- /dev/null +++ b/highs/ipm/hipo/ipm/KktMatrix.cpp @@ -0,0 +1,254 @@ +#include "KktMatrix.h" + +#include "ipm/hipo/auxiliary/Auxiliary.h" + +namespace hipo { + +KktMatrix::KktMatrix(const Model& m, const Regularisation& r, Info& i, + const LogHighs& l) + : model{m}, regul{r}, info{i}, log{l} {} + +Int KktMatrix::buildASstructure() { + // Build lower triangular structure of the augmented system. + // Build values of AS that will not change during the iterations. + + const HighsSparseMatrix& A = model.A(); + const HighsHessian& Q = model.Q(); + const Int n = A.num_col_; + const Int m = A.num_row_; + const Int nzA = A.numNz(); + const Int nzQ = Q.numNz(); + + log.printDevInfo("Building AS structure\n"); + + const Int nzBlock11 = model.qp() ? nzQ : n; + + // AS matrix must fit into HighsInt + if ((Int64)nzBlock11 + m + nzA >= kHighsIInf) return kStatusOverflow; + + ptrAS.resize(n + m + 1); + rowsAS.resize(nzBlock11 + nzA + m); + valAS.resize(nzBlock11 + nzA + m); + + Int next = 0; + + for (Int i = 0; i < n; ++i) { + // diagonal element + rowsAS[next] = i; + next++; + + // column of Q + if (model.qp()) { + assert(Q.index_[Q.start_[i]] == i); + for (Int el = Q.start_[i] + 1; el < Q.start_[i + 1]; ++el) { + rowsAS[next] = Q.index_[el]; + valAS[next] = -Q.value_[el]; // values of AS that will not change + ++next; + } + } + + // column of A + for (Int el = A.start_[i]; el < A.start_[i + 1]; ++el) { + rowsAS[next] = n + A.index_[el]; + valAS[next] = A.value_[el]; // values of AS that will not change + ++next; + } + + ptrAS[i + 1] = next; + } + + // 2,2 block + for (Int i = 0; i < m; ++i) { + rowsAS[next] = n + i; + ++next; + ptrAS[n + i + 1] = ptrAS[n + i] + 1; + } + + return kStatusOk; +} + +Int KktMatrix::buildASvalues(const std::vector& scaling) { + // build AS values that change during iterations. + + assert(!ptrAS.empty() && !rowsAS.empty()); + + const HighsHessian& Q = model.Q(); + const Int n = model.A().num_col_; + + for (Int i = 0; i < n; ++i) { + valAS[ptrAS[i]] = scaling.empty() ? -1.0 : -scaling[i]; + if (model.qp()) valAS[ptrAS[i]] -= model.sense() * model.Q().diag(i); + } + + return kStatusOk; +} + +Int KktMatrix::buildNEstructure() { + // Build lower triangular structure of AAt. + // This approach uses a column-wise copy of A, a partial row-wise copy and a + // vector of corresponding indices. + // NB: A must have sorted columns for this to work + + const HighsSparseMatrix& A = model.A(); + const Int n = A.num_col_; + const Int m = A.num_row_; + const Int nzA = A.numNz(); + + log.printDevInfo("Building NE structure\n"); + + // create partial row-wise representation without values, and array or + // corresponding indices between cw and rw representation + + ptrA_rw.assign(m + 1, 0); + idxA_rw.assign(nzA, 0); + + // pointers of row-start + for (Int el = 0; el < nzA; ++el) ptrA_rw[A.index_[el] + 1]++; + for (Int i = 0; i < m; ++i) ptrA_rw[i + 1] += ptrA_rw[i]; + + std::vector temp = ptrA_rw; + corr_A.assign(nzA, 0); + + // rw-indices and corresponding indices created together + for (Int col = 0; col < n; ++col) { + for (Int el = A.start_[col]; el < A.start_[col + 1]; ++el) { + Int row = A.index_[el]; + + corr_A[temp[row]] = el; + idxA_rw[temp[row]] = col; + temp[row]++; + } + } + + ptrNE.clear(); + rowsNE.clear(); + + // ptr is allocated its exact size + ptrNE.resize(m + 1, 0); + + // keep track if given entry is nonzero, in column considered + std::vector is_nz(m, false); + + // temporary storage of indices + std::vector temp_index(m); + + for (Int row = 0; row < m; ++row) { + // go along the entries of the row, and then down each column. + // this builds the lower triangular part of the row-th column of AAt. + + Int nz_in_col = 0; + + for (Int el = ptrA_rw[row]; el < ptrA_rw[row + 1]; ++el) { + Int col = idxA_rw[el]; + Int corr = corr_A[el]; + + // for each nonzero in the row, go down corresponding column, starting + // from current position + for (Int colEl = corr; colEl < A.start_[col + 1]; ++colEl) { + Int row2 = A.index_[colEl]; + + // row2 is guaranteed to be larger or equal than row + // (provided that the columns of A are sorted) + + // save information that there is nonzero in position (row2,row). + if (!is_nz[row2]) { + is_nz[row2] = true; + temp_index[nz_in_col] = row2; + ++nz_in_col; + } + } + } + // intersection of row with rows below finished. + + // if the total number of nonzeros exceeds the maximum, return error. + if ((Int64)ptrNE[row] + nz_in_col >= + NE_nz_limit.load(std::memory_order_relaxed)) + return kStatusOverflow; + + // update pointers + ptrNE[row + 1] = ptrNE[row] + nz_in_col; + + // now assign indices + for (Int i = 0; i < nz_in_col; ++i) { + Int index = temp_index[i]; + // push_back is better then reserve, because the final length is not known + rowsNE.push_back(index); + is_nz[index] = false; + } + } + + return kStatusOk; +} + +Int KktMatrix::buildNEvalues(const std::vector& scaling) { + // given the NE structure already computed, fill in the NE values + + assert(!ptrNE.empty() && !rowsNE.empty()); + + const HighsSparseMatrix& A = model.A(); + const HighsHessian& Q = model.Q(); + const Int m = A.num_row_; + + valNE.resize(rowsNE.size()); + + std::vector work(m, 0.0); + + for (Int row = 0; row < m; ++row) { + // go along the entries of the row, and then down each column. + // this builds the lower triangular part of the row-th column of AAt. + + for (Int el = ptrA_rw[row]; el < ptrA_rw[row + 1]; ++el) { + Int col = idxA_rw[el]; + Int corr = corr_A[el]; + + double denom = scaling.empty() ? 1.0 : scaling[col]; + denom += regul.primal; + if (model.qp()) denom += model.sense() * Q.diag(col); + + const double mult = 1.0 / denom; + const double row_value = mult * A.value_[corr]; + + // for each nonzero in the row, go down corresponding column, starting + // from current position + for (Int colEl = corr; colEl < A.start_[col + 1]; ++colEl) { + Int row2 = A.index_[colEl]; + + // row2 is guaranteed to be larger or equal than row + // (provided that the columns of A are sorted) + + // compute and accumulate value + double value = row_value * A.value_[colEl]; + work[row2] += value; + } + } + // intersection of row with rows below finished. + + // read from work, using indices of column "row" of AAt + for (Int el = ptrNE[row]; el < ptrNE[row + 1]; ++el) { + Int index = rowsNE[el]; + valNE[el] = work[index]; + work[index] = 0.0; + } + } + + return kStatusOk; +} + +void KktMatrix::freeASmemory() { + // Give up memory used for AS. + freeVector(ptrAS); + freeVector(rowsAS); + freeVector(valAS); +} + +void KktMatrix::freeNEmemory() { + // Give up memory used for NE. + freeVector(ptrNE); + freeVector(rowsNE); + freeVector(valNE); + freeVector(ptrA_rw); + freeVector(idxA_rw); + freeVector(corr_A); +} + +} // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/KktMatrix.h b/highs/ipm/hipo/ipm/KktMatrix.h new file mode 100644 index 0000000000..6d64c89dbd --- /dev/null +++ b/highs/ipm/hipo/ipm/KktMatrix.h @@ -0,0 +1,46 @@ +#ifndef HIPO_KKT_MATRIX_H +#define HIPO_KKT_MATRIX_H + +#include +#include + +#include "Info.h" +#include "Model.h" +#include "ipm/hipo/auxiliary/IntConfig.h" + +namespace hipo { + +struct KktMatrix { + std::vector ptrNE; + std::vector rowsNE; + std::vector valNE; + std::vector ptrA_rw, idxA_rw; + std::vector corr_A; + std::atomic NE_nz_limit{kHighsIInf}; + + std::vector ptrAS; + std::vector rowsAS; + std::vector valAS; + + std::vector iperm; + + const Model& model; + const Regularisation& regul; + Info& info; + const LogHighs& log; + + KktMatrix(const Model& model, const Regularisation& regul, Info& info, + const LogHighs& log); + + Int buildASstructure(); + Int buildASvalues(const std::vector& scaling); + void freeASmemory(); + + Int buildNEstructure(); + Int buildNEvalues(const std::vector& scaling); + void freeNEmemory(); +}; + +} // namespace hipo + +#endif \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/LinearSolver.h b/highs/ipm/hipo/ipm/LinearSolver.h index a44373fff4..f96d5fec6b 100644 --- a/highs/ipm/hipo/ipm/LinearSolver.h +++ b/highs/ipm/hipo/ipm/LinearSolver.h @@ -76,18 +76,6 @@ class LinearSolver { virtual void getReg(std::vector& reg){}; }; -struct KktMatrix { - std::vector ptrNE; - std::vector rowsNE; - std::vector valNE; - - std::vector ptrAS; - std::vector rowsAS; - std::vector valAS; - - std::vector iperm; -}; - } // namespace hipo #endif \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index 67f5b72c10..83f31b5763 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -102,6 +102,10 @@ void Solver::solve() { // iterate object needs to be initialised before potentially interrupting it_.reset(new Iterate(model_, regul_)); + if (!it_) { + info_.status = kStatusError; + return; + } if (checkInterrupt()) return; @@ -132,8 +136,14 @@ bool Solver::initialise() { start_time_ = control_.elapsed(); + kkt_.reset(new KktMatrix(model_, regul_, info_, logH_)); + if (!kkt_) { + info_.status = kStatusError; + return true; + } + // initialise linear solver - LS_.reset(new FactorHiGHSSolver(kkt_, options_, model_, regul_, info_, + LS_.reset(new FactorHiGHSSolver(*kkt_, options_, model_, regul_, info_, it_->data, logH_)); if (Int status = LS_->setup()) { info_.status = (Status)status; diff --git a/highs/ipm/hipo/ipm/Solver.h b/highs/ipm/hipo/ipm/Solver.h index 0fcbd7ec13..a59789610f 100644 --- a/highs/ipm/hipo/ipm/Solver.h +++ b/highs/ipm/hipo/ipm/Solver.h @@ -32,7 +32,7 @@ class Solver { // Linear solver interface std::unique_ptr LS_; - KktMatrix kkt_; + std::unique_ptr kkt_; // Iterate object interface std::unique_ptr it_; From 9f8b67af05a573e973a11b68708366bd2fa77aca Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Wed, 11 Mar 2026 12:23:40 +0000 Subject: [PATCH 04/23] Fix regularisation in kernel --- highs/ipm/hipo/factorhighs/DenseFactKernel.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp index 9a8936421a..d15c950dd0 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp @@ -45,9 +45,9 @@ static void staticReg(double& pivot, Int sign, const Regul& regval, double old_pivot = pivot; if (sign > 0) - pivot += regval.primal; + pivot += regval.dual; else - pivot -= regval.dual; + pivot -= regval.primal; totalreg = pivot - old_pivot; } @@ -111,8 +111,8 @@ static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, assert(r >= 0); res = maxInCol(j, n, r, A, lda); double gamma_r = res.second; - double Arr = sign[r] > 0 ? A[r + lda * r] + regval.primal - : A[r + lda * r] - regval.dual; + double Arr = sign[r] > 0 ? A[r + lda * r] + regval.dual + : A[r + lda * r] - regval.primal; if ((std::abs(Ajj) >= kAlphaBK * gamma_j || std::abs(Ajj) * gamma_r >= kAlphaBK * gamma_j * gamma_j)) { @@ -199,7 +199,7 @@ Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, if (std::isnan(Ajj)) return kRetInvalidPivot; // add regularisation - staticReg(Ajj, pivot_sign[j], regval, totalreg[j]); + //staticReg(Ajj, pivot_sign[j], regval, totalreg[j]); #ifndef HIPO_PIVOTING // add static regularisation From 0e943a463494d19524e2175c0f1c124968963a3a Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 16 Mar 2026 15:28:57 +0000 Subject: [PATCH 05/23] Compute ne nz bounds --- highs/ipm/hipo/ipm/Model.cpp | 28 ++++++++++++++++++++++++++++ highs/ipm/hipo/ipm/Model.h | 7 +++++++ 2 files changed, 35 insertions(+) diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index c1a6fbf6d1..f3d92b606c 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -28,6 +28,7 @@ Int Model::init(const HighsLp& lp, const HighsHessian& Q) { preprocess(); denseColumns(); + nzBounds(); computeNorms(); // double transpose to sort indices of each column @@ -41,6 +42,33 @@ Int Model::init(const HighsLp& lp, const HighsHessian& Q) { return 0; } +void Model::nzBounds() { + // compute lower and upper bounds for the number of nonzeros in normal + // equations. + std::vector mark(m_, false); + NE_nz_lb_ = A_.num_row_; + NE_nz_ub_ = A_.num_row_; + for (Int col = 0; col < A_.num_col_; ++col) { + Int used = 0; + Int unused = 0; + for (Int el = A_.start_[col]; el < A_.start_[col + 1]; ++el) { + const Int row = A_.index_[el]; + if (mark[row]) + ++used; + else { + mark[row] = true; + ++unused; + } + } + + NE_nz_ub_ += (used + unused) * (used + unused - 1) / 2; + NE_nz_lb_ += unused * (unused - 1) / 2 + used * unused; + } + NE_nz_ub_ = std::min(NE_nz_ub_, (Int64)A_.num_row_ * (A_.num_row_ + 1) / 2); + + AS_nz_ = A_.numNz() + A_.num_row_ + (qp() ? Q_.numNz() : A_.num_col_); +} + Int Model::checkData() const { // Check if model provided by the user is ok. // Return kStatusBadModel if something is wrong. diff --git a/highs/ipm/hipo/ipm/Model.h b/highs/ipm/hipo/ipm/Model.h index 11dec6e4bd..c1eace849f 100644 --- a/highs/ipm/hipo/ipm/Model.h +++ b/highs/ipm/hipo/ipm/Model.h @@ -49,8 +49,10 @@ class Model { HighsSparseMatrix A_{}; HighsHessian Q_{}; std::vector constraints_{}; + Int num_dense_cols_{}; double max_col_density_{}; + Int64 AS_nz_, NE_nz_lb_, NE_nz_ub_; std::vector colscale_, rowscale_; @@ -89,6 +91,8 @@ class Model { double normUnscaledObj() const { return norm_unscaled_obj_; } double normUnscaledRhs() const { return norm_unscaled_rhs_; } + void nzBounds(); + // Check if variable has finite lower/upper bound bool hasLb(Int j) const { return std::isfinite(lower_[j]); } bool hasUb(Int j) const { return std::isfinite(upper_[j]); } @@ -121,6 +125,9 @@ class Model { double oneNormCols(Int i) const { return one_norm_cols_[i]; } double infNormRows(Int i) const { return inf_norm_rows_[i]; } double infNormCols(Int i) const { return inf_norm_cols_[i]; } + Int64 nzAS() const { return AS_nz_; } + Int64 nzNElb() const { return NE_nz_lb_; } + Int64 nzNEub() const { return NE_nz_ub_; } Int loadIntoIpx(ipx::LpSolver& lps) const; From beb267a4cc5e4d6b6e92cff81056c41610ea447a Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 31 Mar 2026 13:43:00 +0100 Subject: [PATCH 06/23] Fix timers --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 19 ++----------------- highs/ipm/hipo/ipm/KktMatrix.cpp | 11 +++++++++++ 2 files changed, 13 insertions(+), 17 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 4562a2baf3..b7272798c0 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -86,17 +86,12 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { // only execute factorisation if it has not been done yet assert(!this->valid_); - Clock clock; - - // build matrix kkt_.buildASvalues(scaling); - info_.matrix_time += clock.stop(); // set static regularisation, since it may have changed FH_.setRegularisation(regul_.primal, regul_.dual); - // factorise matrix - clock.start(); + Clock clock; if (FH_.factorise(S_, kkt_.rowsAS, kkt_.ptrAS, kkt_.valAS)) return kStatusErrorFactorise; info_.factor_time += clock.stop(); @@ -110,20 +105,14 @@ Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { // only execute factorisation if it has not been done yet assert(!this->valid_); - Clock clock; - - // build matrix kkt_.buildNEvalues(scaling); - info_.matrix_time += clock.stop(); // set static regularisation, since it may have changed FH_.setRegularisation(regul_.primal, regul_.dual); - // factorise - clock.start(); + Clock clock; if (FH_.factorise(S_, kkt_.rowsNE, kkt_.ptrNE, kkt_.valNE)) return kStatusErrorFactorise; - info_.factor_time += clock.stop(); info_.factor_number++; @@ -219,9 +208,7 @@ Int FactorHiGHSSolver::chooseNla() { model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; } else { - Clock clock; Int status = kkt_.buildNEstructure(); - info_.NE_structure_time = clock.stop(); if (status) { failure_NE = true; if (status == kStatusOverflow) { @@ -338,8 +325,6 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // - If ordering is "amd", "metis", "rcm" run only the ordering requested. // - If ordering is "choose", run "amd", "metis", and choose the best. - Clock clock; - // select which fill-reducing orderings should be tried std::vector orderings_to_try; if (options_.ordering != kHighsChooseString) diff --git a/highs/ipm/hipo/ipm/KktMatrix.cpp b/highs/ipm/hipo/ipm/KktMatrix.cpp index 8dc050ed56..1926c4c337 100644 --- a/highs/ipm/hipo/ipm/KktMatrix.cpp +++ b/highs/ipm/hipo/ipm/KktMatrix.cpp @@ -11,6 +11,7 @@ KktMatrix::KktMatrix(const Model& m, const Regularisation& r, Info& i, Int KktMatrix::buildASstructure() { // Build lower triangular structure of the augmented system. // Build values of AS that will not change during the iterations. + Clock clock; const HighsSparseMatrix& A = model.A(); const HighsHessian& Q = model.Q(); @@ -64,6 +65,8 @@ Int KktMatrix::buildASstructure() { ptrAS[n + i + 1] = ptrAS[n + i] + 1; } + info.AS_structure_time = clock.stop(); + return kStatusOk; } @@ -71,6 +74,7 @@ Int KktMatrix::buildASvalues(const std::vector& scaling) { // build AS values that change during iterations. assert(!ptrAS.empty() && !rowsAS.empty()); + Clock clock; const HighsHessian& Q = model.Q(); const Int n = model.A().num_col_; @@ -80,6 +84,8 @@ Int KktMatrix::buildASvalues(const std::vector& scaling) { if (model.qp()) valAS[ptrAS[i]] -= model.sense() * model.Q().diag(i); } + info.matrix_time += clock.stop(); + return kStatusOk; } @@ -88,6 +94,7 @@ Int KktMatrix::buildNEstructure() { // This approach uses a column-wise copy of A, a partial row-wise copy and a // vector of corresponding indices. // NB: A must have sorted columns for this to work + Clock clock; const HighsSparseMatrix& A = model.A(); const Int n = A.num_col_; @@ -177,6 +184,7 @@ Int KktMatrix::buildNEstructure() { } } + info.NE_structure_time = clock.stop(); return kStatusOk; } @@ -184,6 +192,7 @@ Int KktMatrix::buildNEvalues(const std::vector& scaling) { // given the NE structure already computed, fill in the NE values assert(!ptrNE.empty() && !rowsNE.empty()); + Clock clock; const HighsSparseMatrix& A = model.A(); const HighsHessian& Q = model.Q(); @@ -231,6 +240,8 @@ Int KktMatrix::buildNEvalues(const std::vector& scaling) { } } + info.matrix_time += clock.stop(); + return kStatusOk; } From 2fcc9228261e634faa3ac9a6e7c38fa2088aef44 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 16 Mar 2026 16:04:47 +0000 Subject: [PATCH 07/23] Use bounds to predict as or ne failure --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 41 +++++++++++++++--------- highs/ipm/hipo/ipm/Model.cpp | 9 +++++- highs/ipm/hipo/ipm/Parameters.h | 1 + 3 files changed, 35 insertions(+), 16 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index b7272798c0..2f7597ad38 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -203,9 +203,13 @@ Int FactorHiGHSSolver::chooseNla() { bool overflow_AS = false; auto run_structure_NE = [&]() { - if ((model_.m() > kMinRowsForDensity && - model_.maxColDensity() > kDenseColThresh) || - model_.nonSeparableQp() || model_.m() == 0) { + bool has_dense_cols = model_.m() > kMinRowsForDensity && + model_.maxColDensity() > kDenseColThresh; + bool expect_AS_much_cheaper = + model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; + + if (has_dense_cols || expect_AS_much_cheaper || model_.nonSeparableQp() || + model_.m() == 0) { failure_NE = true; } else { Int status = kkt_.buildNEstructure(); @@ -227,19 +231,26 @@ Int FactorHiGHSSolver::chooseNla() { }; auto run_analyse_AS = [&]() { - Int AS_status = kkt_.buildASstructure(); - if (!AS_status) AS_status = analyseAS(symb_AS); - if (AS_status) failure_AS = true; - if (AS_status == kStatusOverflow) { - log_.printDevInfo("Integer overflow forming AS matrix\n"); - overflow_AS = true; - } + bool expect_NE_much_cheaper = + model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; + + if (expect_NE_much_cheaper) + failure_AS = true; + else { + Int AS_status = kkt_.buildASstructure(); + if (!AS_status) AS_status = analyseAS(symb_AS); + if (AS_status) failure_AS = true; + if (AS_status == kStatusOverflow) { + log_.printDevInfo("Integer overflow forming AS matrix\n"); + overflow_AS = true; + } - // If NE has more nonzeros than the factor of AS, then it's likely that AS - // will be preferred, so stop computation of NE. - Int64 NE_nz_limit = symb_AS.nz() * kSymbNzMult; - if (failure_AS || NE_nz_limit > kHighsIInf) NE_nz_limit = kHighsIInf; - kkt_.NE_nz_limit.store(NE_nz_limit, std::memory_order_relaxed); + // If NE has more nonzeros than the factor of AS, then it's likely that AS + // will be preferred, so stop computation of NE. + Int64 NE_nz_limit = symb_AS.nz() * kSymbNzMult; + if (failure_AS || NE_nz_limit > kHighsIInf) NE_nz_limit = kHighsIInf; + kkt_.NE_nz_limit.store(NE_nz_limit, std::memory_order_relaxed); + } }; // In parallel, run AS analyse and build NE structure. NE analyse runs only diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index f3d92b606c..406fca1b17 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -343,7 +343,14 @@ void Model::print(const LogHighs& log) const { log_stream << textline("Scaling coefficients:") << "[" << sci(scalemin, 5, 1) << ", " << sci(scalemax, 5, 1) << "]\n"; - if (log.debug(1)) preprocessor_.print(log_stream); + if (log.debug(1)) { + preprocessor_.print(log_stream); + + log_stream << "Expected nnz: "; + log_stream << "AS " << sci(AS_nz_, 0, 1) << "; "; + log_stream << "NE [" << sci(NE_nz_lb_, 0, 1) << ", " << sci(NE_nz_ub_, 0, 1) + << "]\n"; + } log.print(log_stream); } diff --git a/highs/ipm/hipo/ipm/Parameters.h b/highs/ipm/hipo/ipm/Parameters.h index af2c504d9c..4c6a07c760 100644 --- a/highs/ipm/hipo/ipm/Parameters.h +++ b/highs/ipm/hipo/ipm/Parameters.h @@ -42,6 +42,7 @@ const double kFlopsOrderingThresh = 1.2; // parameters for dense columns const double kDenseColThresh = 0.5; const Int kMinRowsForDensity = 2000; +const double kNzBoundsRatio = 50.0; // parameters for iterative refinement const Int kMaxIterRefine = 3; From 94ffd99781af8f86ccff49c7f9f94fb4a7129f75 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Thu, 19 Mar 2026 14:37:36 +0000 Subject: [PATCH 08/23] Print free vars --- highs/ipm/hipo/ipm/Model.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index 406fca1b17..96f2646cb2 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -280,12 +280,15 @@ void Model::print(const LogHighs& log) const { // compute max and min for bounds intervals double boundmin = kHighsInf; double boundmax = 0.0; + Int free_vars = 0; for (Int i = 0; i < n_; ++i) { if (std::isfinite(lower_[i]) && std::isfinite(upper_[i])) { const double diff = std::abs(upper_[i] - lower_[i]); boundmin = std::min(boundmin, diff); boundmax = std::max(boundmax, diff); } + + if (!std::isfinite(lower_[i]) && !std::isfinite(upper_[i])) free_vars++; } if (std::isinf(boundmin)) boundmin = 0.0; @@ -345,6 +348,7 @@ void Model::print(const LogHighs& log) const { if (log.debug(1)) { preprocessor_.print(log_stream); + log_stream << textline("Free variables:") << integer(free_vars) << '\n'; log_stream << "Expected nnz: "; log_stream << "AS " << sci(AS_nz_, 0, 1) << "; "; From ed18615d5422d236ca649e6357ccd495548fbd07 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 31 Mar 2026 14:18:02 +0100 Subject: [PATCH 09/23] Print message when AS or NE is skipped due to nz bounds --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 2f7597ad38..22f83f12c5 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -211,6 +211,7 @@ Int FactorHiGHSSolver::chooseNla() { if (has_dense_cols || expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; + log_.printDevInfo("NE skipped due to nz bounds\n"); } else { Int status = kkt_.buildNEstructure(); if (status) { @@ -234,9 +235,10 @@ Int FactorHiGHSSolver::chooseNla() { bool expect_NE_much_cheaper = model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; - if (expect_NE_much_cheaper) + if (expect_NE_much_cheaper) { failure_AS = true; - else { + log_.printDevInfo("AS skipped due to nz bounds\n"); + } else { Int AS_status = kkt_.buildASstructure(); if (!AS_status) AS_status = analyseAS(symb_AS); if (AS_status) failure_AS = true; From aa48ee2b575f8d899ba889135f0f75021c457d73 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 31 Mar 2026 14:46:59 +0100 Subject: [PATCH 10/23] Do not skip AS if NE not possible --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 22f83f12c5..41cbef244a 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -211,7 +211,7 @@ Int FactorHiGHSSolver::chooseNla() { if (has_dense_cols || expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; - log_.printDevInfo("NE skipped due to nz bounds\n"); + log_.printDevInfo("NE skipped\n"); } else { Int status = kkt_.buildNEstructure(); if (status) { @@ -235,9 +235,11 @@ Int FactorHiGHSSolver::chooseNla() { bool expect_NE_much_cheaper = model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; - if (expect_NE_much_cheaper) { + bool can_skip_AS = !(model_.nonSeparableQp() || model_.m() == 0); + + if (expect_NE_much_cheaper && can_skip_AS) { failure_AS = true; - log_.printDevInfo("AS skipped due to nz bounds\n"); + log_.printDevInfo("AS skipped\n"); } else { Int AS_status = kkt_.buildASstructure(); if (!AS_status) AS_status = analyseAS(symb_AS); From db9a02567f5918dd549acb355eaad95e67516aed Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 31 Mar 2026 14:51:24 +0100 Subject: [PATCH 11/23] Remove dense columns stuff --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 4 +--- highs/ipm/hipo/ipm/Info.h | 4 ---- highs/ipm/hipo/ipm/Model.cpp | 20 +------------------- highs/ipm/hipo/ipm/Model.h | 5 ----- highs/ipm/hipo/ipm/Parameters.h | 4 +--- highs/ipm/hipo/ipm/Solver.cpp | 3 --- 6 files changed, 3 insertions(+), 37 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 41cbef244a..a2653a9d80 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -203,12 +203,10 @@ Int FactorHiGHSSolver::chooseNla() { bool overflow_AS = false; auto run_structure_NE = [&]() { - bool has_dense_cols = model_.m() > kMinRowsForDensity && - model_.maxColDensity() > kDenseColThresh; bool expect_AS_much_cheaper = model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; - if (has_dense_cols || expect_AS_much_cheaper || model_.nonSeparableQp() || + if ( expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; log_.printDevInfo("NE skipped\n"); diff --git a/highs/ipm/hipo/ipm/Info.h b/highs/ipm/hipo/ipm/Info.h index 0f17677af1..14ef24138d 100644 --- a/highs/ipm/hipo/ipm/Info.h +++ b/highs/ipm/hipo/ipm/Info.h @@ -47,10 +47,6 @@ struct Info { // Counters Int factor_number{}; Int solve_number{}; - - // Information on dense columns - Int num_dense_cols{}; - double max_col_density{}; }; } // namespace hipo diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index 96f2646cb2..e61acf3655 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -27,7 +27,6 @@ Int Model::init(const HighsLp& lp, const HighsHessian& Q) { preprocess(); - denseColumns(); nzBounds(); computeNorms(); @@ -184,9 +183,7 @@ void Model::print(const LogHighs& log) const { log_stream << textline("Rows:") << sci(m_, 0, 1) << '\n'; log_stream << textline("Cols:") << sci(n_, 0, 1) << '\n'; log_stream << textline("Nnz A:") << sci(A_.numNz(), 0, 1) << '\n'; - if (num_dense_cols_ > 0) - log_stream << textline("Dense cols:") << integer(num_dense_cols_, 0) - << '\n'; + if (qp()) { log_stream << textline("Nnz Q:") << sci(Q_.numNz(), 0, 1); if (nonSeparableQp()) @@ -359,21 +356,6 @@ void Model::print(const LogHighs& log) const { log.print(log_stream); } -void Model::denseColumns() { - // Compute the maximum density of any column of A and count the number of - // dense columns. - - max_col_density_ = 0.0; - num_dense_cols_ = 0; - for (Int col = 0; col < n_; ++col) { - Int col_nz = A_.start_[col + 1] - A_.start_[col]; - double col_density = (double)col_nz / m_; - max_col_density_ = std::max(max_col_density_, col_density); - if (A_.num_row_ > kMinRowsForDensity && col_density > kDenseColThresh) - ++num_dense_cols_; - } -} - Int Model::loadIntoIpx(ipx::LpSolver& lps) const { Int ipx_m, ipx_n; std::vector ipx_b, ipx_c, ipx_lower, ipx_upper, ipx_A_vals; diff --git a/highs/ipm/hipo/ipm/Model.h b/highs/ipm/hipo/ipm/Model.h index c1eace849f..e36541e50c 100644 --- a/highs/ipm/hipo/ipm/Model.h +++ b/highs/ipm/hipo/ipm/Model.h @@ -50,8 +50,6 @@ class Model { HighsHessian Q_{}; std::vector constraints_{}; - Int num_dense_cols_{}; - double max_col_density_{}; Int64 AS_nz_, NE_nz_lb_, NE_nz_ub_; std::vector colscale_, rowscale_; @@ -68,7 +66,6 @@ class Model { inf_norm_rows_; void preprocess(); - void denseColumns(); Int checkData() const; void computeNorms(); @@ -119,8 +116,6 @@ class Model { bool ready() const { return ready_; } bool scaled() const { return !colscale_.empty(); } double offset() const { return offset_; } - double maxColDensity() const { return max_col_density_; } - Int numDenseCols() const { return num_dense_cols_; } double oneNormRows(Int i) const { return one_norm_rows_[i]; } double oneNormCols(Int i) const { return one_norm_cols_[i]; } double infNormRows(Int i) const { return inf_norm_rows_[i]; } diff --git a/highs/ipm/hipo/ipm/Parameters.h b/highs/ipm/hipo/ipm/Parameters.h index 4c6a07c760..c80f61b230 100644 --- a/highs/ipm/hipo/ipm/Parameters.h +++ b/highs/ipm/hipo/ipm/Parameters.h @@ -39,9 +39,7 @@ const double kMaxTreeDepth = 1000; // parameters for choice of ordering const double kFlopsOrderingThresh = 1.2; -// parameters for dense columns -const double kDenseColThresh = 0.5; -const Int kMinRowsForDensity = 2000; +// parameters for skipping AS or NE const double kNzBoundsRatio = 50.0; // parameters for iterative refinement diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index 83f31b5763..4334825ba6 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -174,9 +174,6 @@ bool Solver::initialise() { void Solver::terminate() { info_.ipm_iter = iter_; if (info_.status == kStatusNotRun) info_.status = kStatusMaxIter; - - info_.num_dense_cols = model_.numDenseCols(); - info_.max_col_density = model_.maxColDensity(); } bool Solver::prepareIter() { From 5e4e84213237de58ac621e2b868e6a780573bd2f Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 31 Mar 2026 16:42:14 +0100 Subject: [PATCH 12/23] Fix bad infeasibility detection --- highs/ipm/hipo/ipm/Iterate.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/highs/ipm/hipo/ipm/Iterate.cpp b/highs/ipm/hipo/ipm/Iterate.cpp index a2ca24ea33..84cbd35579 100644 --- a/highs/ipm/hipo/ipm/Iterate.cpp +++ b/highs/ipm/hipo/ipm/Iterate.cpp @@ -619,12 +619,8 @@ bool Iterate::stagnation(std::stringstream& log_stream) { // infeasibilities jumping back up const double thresh_inf_to_best = 1e12; - best_pinf_ = std::min(best_pinf_, pinf); - best_dinf_ = std::min(best_dinf_, dinf); - - // if the best is zero, the test would always be triggered - best_pinf_ = std::max(best_pinf_, 1e-16); - best_dinf_ = std::max(best_dinf_, 1e-16); + if (pinf > 0.0) best_pinf_ = std::min(best_pinf_, pinf); + if (dinf > 0.0) best_dinf_ = std::min(best_dinf_, dinf); if (pinf > thresh_inf_to_best * best_pinf_ || dinf > thresh_inf_to_best * best_dinf_) { From 67a4d90e607b51cadb2d51beee0500816dd82add Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 6 Apr 2026 11:56:19 +0100 Subject: [PATCH 13/23] Print analyse time --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index a2653a9d80..61d946a722 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -177,9 +177,15 @@ Int FactorHiGHSSolver::solveNE(const std::vector& rhs, // ========================================================================= Int FactorHiGHSSolver::setup() { + Clock clock; + if (Int status = setNla()) return status; setParallel(); + std::stringstream log_stream; + log_stream << textline("Analyse time:") << fix(clock.stop(), 0, 2) << '\n'; + log_.print(log_stream); + S_.print(log_, log_.debug(1)); // Warn about large memory consumption @@ -206,8 +212,7 @@ Int FactorHiGHSSolver::chooseNla() { bool expect_AS_much_cheaper = model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; - if ( expect_AS_much_cheaper || model_.nonSeparableQp() || - model_.m() == 0) { + if (expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; log_.printDevInfo("NE skipped\n"); } else { @@ -343,8 +348,8 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, if (options_.ordering != kHighsChooseString) orderings_to_try.push_back(options_.ordering); else { - orderings_to_try.push_back("amd"); - orderings_to_try.push_back("metis"); + orderings_to_try.push_back(kHipoAmdString); + orderings_to_try.push_back(kHipoMetisString); // rcm is much worse in general, so no point in trying for now } From 3fe4bcb92ab01d6df4f934b5756f1a562d311146 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Fri, 10 Apr 2026 11:24:55 +0100 Subject: [PATCH 14/23] Simplify logging --- cmake/sources.cmake | 6 +- highs/ipm/hipo/auxiliary/KrylovMethods.cpp | 2 - highs/ipm/hipo/auxiliary/Log.cpp | 60 ------- highs/ipm/hipo/auxiliary/Log.h | 62 -------- highs/ipm/hipo/auxiliary/Logger.cpp | 44 ++++++ highs/ipm/hipo/auxiliary/Logger.h | 80 ++++++++++ highs/ipm/hipo/factorhighs/Analyse.cpp | 10 +- highs/ipm/hipo/factorhighs/Analyse.h | 6 +- highs/ipm/hipo/factorhighs/DataCollector.cpp | 10 +- highs/ipm/hipo/factorhighs/DataCollector.h | 6 +- .../ipm/hipo/factorhighs/DenseFactHybrid.cpp | 1 - .../ipm/hipo/factorhighs/DenseFactKernel.cpp | 1 - highs/ipm/hipo/factorhighs/FactorHiGHS.cpp | 24 +-- highs/ipm/hipo/factorhighs/FactorHiGHS.h | 14 +- highs/ipm/hipo/factorhighs/Factorise.cpp | 24 +-- highs/ipm/hipo/factorhighs/Factorise.h | 6 +- highs/ipm/hipo/factorhighs/Numeric.cpp | 1 - highs/ipm/hipo/factorhighs/Symbolic.cpp | 5 +- highs/ipm/hipo/factorhighs/Symbolic.h | 4 +- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 76 ++++----- highs/ipm/hipo/ipm/FactorHiGHSSolver.h | 4 +- highs/ipm/hipo/ipm/KktMatrix.cpp | 8 +- highs/ipm/hipo/ipm/KktMatrix.h | 4 +- highs/ipm/hipo/ipm/LogHighs.cpp | 71 --------- highs/ipm/hipo/ipm/LogHighs.h | 33 ---- highs/ipm/hipo/ipm/Model.cpp | 10 +- highs/ipm/hipo/ipm/Model.h | 4 +- highs/ipm/hipo/ipm/Solver.cpp | 149 ++++++++---------- highs/ipm/hipo/ipm/Solver.h | 4 +- 29 files changed, 298 insertions(+), 431 deletions(-) delete mode 100644 highs/ipm/hipo/auxiliary/Log.cpp delete mode 100644 highs/ipm/hipo/auxiliary/Log.h create mode 100644 highs/ipm/hipo/auxiliary/Logger.cpp create mode 100644 highs/ipm/hipo/auxiliary/Logger.h delete mode 100644 highs/ipm/hipo/ipm/LogHighs.cpp delete mode 100644 highs/ipm/hipo/ipm/LogHighs.h diff --git a/cmake/sources.cmake b/cmake/sources.cmake index ac872418eb..5074a4802b 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -183,7 +183,6 @@ set(hipo_sources ipm/hipo/ipm/Control.cpp ipm/hipo/ipm/Iterate.cpp ipm/hipo/ipm/KktMatrix.cpp - ipm/hipo/ipm/LogHighs.cpp ipm/hipo/ipm/Model.cpp ipm/hipo/ipm/PreProcess.cpp ipm/hipo/ipm/Refine.cpp @@ -198,7 +197,6 @@ set(hipo_headers ipm/hipo/ipm/Iterate.h ipm/hipo/ipm/KktMatrix.h ipm/hipo/ipm/LinearSolver.h - ipm/hipo/ipm/LogHighs.h ipm/hipo/ipm/Model.h ipm/hipo/ipm/Options.h ipm/hipo/ipm/PreProcess.h @@ -248,14 +246,14 @@ set(factor_highs_headers set(hipo_util_sources ipm/hipo/auxiliary/Auxiliary.cpp ipm/hipo/auxiliary/KrylovMethods.cpp - ipm/hipo/auxiliary/Log.cpp + ipm/hipo/auxiliary/Logger.cpp ipm/hipo/auxiliary/VectorOperations.cpp) set(hipo_util_headers ipm/hipo/auxiliary/Auxiliary.h ipm/hipo/auxiliary/IntConfig.h ipm/hipo/auxiliary/KrylovMethods.h - ipm/hipo/auxiliary/Log.h + ipm/hipo/auxiliary/Logger.h ipm/hipo/auxiliary/mycblas.h ipm/hipo/auxiliary/OrderingPrint.h ipm/hipo/auxiliary/VectorOperations.h) diff --git a/highs/ipm/hipo/auxiliary/KrylovMethods.cpp b/highs/ipm/hipo/auxiliary/KrylovMethods.cpp index b4f79fd75f..5b20e29812 100644 --- a/highs/ipm/hipo/auxiliary/KrylovMethods.cpp +++ b/highs/ipm/hipo/auxiliary/KrylovMethods.cpp @@ -4,7 +4,6 @@ #include #include "VectorOperations.h" -#include "ipm/hipo/auxiliary/Log.h" namespace hipo { @@ -182,7 +181,6 @@ Int Cg(const AbstractMatrix* M, const AbstractMatrix* P, ++iter; if (isNanVector(x)) { - // log->printDevInfo("CG: x is nan at iter %d\n", (int)iter); break; } } diff --git a/highs/ipm/hipo/auxiliary/Log.cpp b/highs/ipm/hipo/auxiliary/Log.cpp deleted file mode 100644 index 742939bf50..0000000000 --- a/highs/ipm/hipo/auxiliary/Log.cpp +++ /dev/null @@ -1,60 +0,0 @@ -#include "Log.h" - -namespace hipo { - -Log::Log(Int level) : dev_level_{level} {} - -void Log::print(std::stringstream& ss) const { printf("%s", ss.str().c_str()); } -void Log::printw(std::stringstream& ss) const { - printf("WARNING: %s", ss.str().c_str()); -} -void Log::printe(std::stringstream& ss) const { - printf("ERROR: %s", ss.str().c_str()); -} -void Log::print(const char* c) const { printf("%s", c); } -void Log::printw(const char* c) const { printf("WARNING: %s", c); } -void Log::printe(const char* c) const { printf("ERROR: %s", c); } -void Log::printDevInfo(std::stringstream& ss) const { - if (dev_level_ >= 1) printf("%s", ss.str().c_str()); -} -void Log::printDevDetailed(std::stringstream& ss) const { - if (dev_level_ >= 2) printf("%s", ss.str().c_str()); -} -void Log::printDevVerbose(std::stringstream& ss) const { - if (dev_level_ >= 3) printf("%s", ss.str().c_str()); -} -void Log::printDevInfo(const char* c) const { - if (dev_level_ >= 1) printf("%s", c); -} -void Log::printDevDetailed(const char* c) const { - if (dev_level_ >= 2) printf("%s", c); -} -void Log::printDevVerbose(const char* c) const { - if (dev_level_ >= 3) printf("%s", c); -} -bool Log::debug(Int level) const { - // Return true if level agrees with dev_level - if (dev_level_ == 1) return level == 1; - if (dev_level_ == 2) return level == 1 || level == 2; - if (dev_level_ == 3) return level == 1 || level == 2 || level == 3; - return false; -} - -std::string format(double d, Int width, Int prec, - std::ios_base::fmtflags floatfield) { - std::ostringstream s; - s.precision(prec); - s.width(width); - s.setf(floatfield, std::ios_base::floatfield); - s << d; - return s.str(); -} - -std::string integer(Int i, Int width) { - std::ostringstream s; - s.width(width); - s << i; - return s.str(); -} - -} // namespace hipo diff --git a/highs/ipm/hipo/auxiliary/Log.h b/highs/ipm/hipo/auxiliary/Log.h deleted file mode 100644 index 72ee056951..0000000000 --- a/highs/ipm/hipo/auxiliary/Log.h +++ /dev/null @@ -1,62 +0,0 @@ -#ifndef HIPO_LOG_H -#define HIPO_LOG_H - -#include - -#include "io/HighsIO.h" -#include "ipm/hipo/auxiliary/IntConfig.h" - -// Base class for logging. -// Use print for normal logging, printw for warnings, printe for errors. -// dev_level = 0 does not print any debug information. -// dev_level = 1 prints debug level Info. -// dev_level = 2 prints debug levels Info, Detailed. -// dev_level = 3 prints debug levels Info, Detailed, Verbose. - -namespace hipo { - -class Log { - Int dev_level_ = 0; - - public: - Log(Int level = 0); - - virtual void print(std::stringstream& ss) const; - virtual void printw(std::stringstream& ss) const; - virtual void printe(std::stringstream& ss) const; - virtual void print(const char* c) const; - virtual void printw(const char* c) const; - virtual void printe(const char* c) const; - - virtual void printDevInfo(std::stringstream& ss) const; - virtual void printDevDetailed(std::stringstream& ss) const; - virtual void printDevVerbose(std::stringstream& ss) const; - virtual void printDevInfo(const char* c) const; - virtual void printDevDetailed(const char* c) const; - virtual void printDevVerbose(const char* c) const; - - virtual bool debug(Int level) const; -}; - -// Functions to print using streams, taken from IPX. -std::string format(double d, Int width, Int prec, - std::ios_base::fmtflags floatfield); -std::string integer(Int i, Int width = 0); -inline std::string sci(double d, Int width, Int prec) { - return format(d, width, prec, std::ios_base::scientific); -} -inline std::string fix(double d, Int width, Int prec) { - return format(d, width, prec, std::ios_base::fixed); -} -template -std::string textline(const T& text) { - std::ostringstream s; - s.setf(std::ios_base::left); - s.width(32); - s << text; - return s.str(); -} - -} // namespace hipo - -#endif \ No newline at end of file diff --git a/highs/ipm/hipo/auxiliary/Logger.cpp b/highs/ipm/hipo/auxiliary/Logger.cpp new file mode 100644 index 0000000000..3cc4688d44 --- /dev/null +++ b/highs/ipm/hipo/auxiliary/Logger.cpp @@ -0,0 +1,44 @@ +#include "Logger.h" + +namespace hipo { + +Logger::Logger(Int level) : dev_level_{level} {} + +void Logger::setOptions(const HighsLogOptions* log_options) { + log_options_ = log_options; +} + +bool Logger::debug(Int level) const { + // Return true if level agrees with log_dev_level + + if (log_options_ && log_options_->log_dev_level) { + if (*log_options_->log_dev_level == kHighsLogDevLevelInfo) + return level == 1; + + if (*log_options_->log_dev_level == kHighsLogDevLevelDetailed) + return level == 1 || level == 2; + + if (*log_options_->log_dev_level == kHighsLogDevLevelVerbose) + return level == 1 || level == 2 || level == 3; + } + return false; +} + +std::string format(double d, Int width, Int prec, + std::ios_base::fmtflags floatfield) { + std::ostringstream s; + s.precision(prec); + s.width(width); + s.setf(floatfield, std::ios_base::floatfield); + s << d; + return s.str(); +} + +std::string integer(Int i, Int width) { + std::ostringstream s; + s.width(width); + s << i; + return s.str(); +} + +} // namespace hipo diff --git a/highs/ipm/hipo/auxiliary/Logger.h b/highs/ipm/hipo/auxiliary/Logger.h new file mode 100644 index 0000000000..b670c90ee0 --- /dev/null +++ b/highs/ipm/hipo/auxiliary/Logger.h @@ -0,0 +1,80 @@ +#ifndef HIPO_LOGGER_H +#define HIPO_LOGGER_H + +#include + +#include "io/HighsIO.h" +#include "ipm/hipo/auxiliary/IntConfig.h" + +// Class for Highs logging +// Use setOptions to pass the log_options. + +namespace hipo { + +class Logger { + const Int dev_level_ = 0; + const HighsLogOptions* log_options_; + + public: + Logger(Int level = 0); + void setOptions(const HighsLogOptions* log_options); + bool debug(Int level) const; + + template + void print(const char* format, Args... args) const { + if (log_options_) + highsLogUser(*log_options_, HighsLogType::kInfo, format, args...); + } + + template + void printw(const char* format, Args... args) const { + if (log_options_) + highsLogUser(*log_options_, HighsLogType::kWarning, format, args...); + } + + template + void printe(const char* format, Args... args) const { + if (log_options_) + highsLogUser(*log_options_, HighsLogType::kError, format, args...); + } + + template + void printInfo(const char* format, Args... args) const { + if (log_options_) + highsLogDev(*log_options_, HighsLogType::kInfo, format, args...); + } + + template + void printDetailed(const char* format, Args... args) const { + if (log_options_) + highsLogDev(*log_options_, HighsLogType::kDetailed, format, args...); + } + + template + void printVerbose(const char* format, Args... args) const { + if (log_options_) + highsLogDev(*log_options_, HighsLogType::kVerbose, format, args...); + } +}; + +// Functions to print using streams, taken from IPX. +std::string format(double d, Int width, Int prec, + std::ios_base::fmtflags floatfield); +std::string integer(Int i, Int width = 0); +inline std::string sci(double d, Int width, Int prec) { + return format(d, width, prec, std::ios_base::scientific); +} +inline std::string fix(double d, Int width, Int prec) { + return format(d, width, prec, std::ios_base::fixed); +} +template +std::string textline(const T& text) { + std::ostringstream s; + s.setf(std::ios_base::left); + s.width(32); + s << text; + return s.str(); +} + +} // namespace hipo +#endif \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/Analyse.cpp b/highs/ipm/hipo/factorhighs/Analyse.cpp index f4153e72a2..73a4bdb593 100644 --- a/highs/ipm/hipo/factorhighs/Analyse.cpp +++ b/highs/ipm/hipo/factorhighs/Analyse.cpp @@ -12,7 +12,7 @@ #include "ReturnValues.h" #include "amd/amd.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "metis/metis.h" #include "rcm/rcm.h" @@ -22,9 +22,9 @@ const Int64 int32_limit = std::numeric_limits::max(); const Int64 int64_limit = std::numeric_limits::max(); Analyse::Analyse(const std::vector& rows, const std::vector& ptr, - const std::vector& signs, Int nb, const Log* log, + const std::vector& signs, Int nb, const Logger* logger, DataCollector& data, const std::vector& perm) - : log_{log}, data_{data} { + : logger_{logger}, data_{data} { // Input the symmetric matrix to be analysed in CSC format. // rows contains the row indices. // ptr contains the starting points of each column. @@ -882,7 +882,7 @@ void Analyse::relativeIndClique() { } else if (consecutive_sums_[sn][i] == 1) { consecutive_sums_[sn][i] = consecutive_sums_[sn][i + 1] + 1; } else { - if (log_) log_->printDevInfo("Error in consecutiveSums\n"); + if (logger_) logger_->printInfo("Error in consecutiveSums\n"); } } } @@ -1282,7 +1282,7 @@ Int Analyse::run(Symbolic& S) { } if (checkOverflow()) { - if (log_) log_->printe("Integer overflow in analyse phase\n"); + if (logger_) logger_->printe("Integer overflow in analyse phase\n"); return kRetIntOverflow; } diff --git a/highs/ipm/hipo/factorhighs/Analyse.h b/highs/ipm/hipo/factorhighs/Analyse.h index c82b4c5f6b..330578615d 100644 --- a/highs/ipm/hipo/factorhighs/Analyse.h +++ b/highs/ipm/hipo/factorhighs/Analyse.h @@ -7,7 +7,7 @@ #include "DataCollector.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { @@ -82,7 +82,7 @@ class Analyse { // block size Int nb_{}; - const Log* log_; + const Logger* logger_; DataCollector& data_; // Functions to perform analyse phase @@ -110,7 +110,7 @@ class Analyse { public: // Constructor: matrix must be in lower triangular format Analyse(const std::vector& rows, const std::vector& ptr, - const std::vector& signs, Int nb, const Log* log, + const std::vector& signs, Int nb, const Logger* logger, DataCollector& data, const std::vector& perm); // Run analyse phase and save the result in Symbolic object S diff --git a/highs/ipm/hipo/factorhighs/DataCollector.cpp b/highs/ipm/hipo/factorhighs/DataCollector.cpp index 62d132a8c5..c0e70d4801 100644 --- a/highs/ipm/hipo/factorhighs/DataCollector.cpp +++ b/highs/ipm/hipo/factorhighs/DataCollector.cpp @@ -1,7 +1,7 @@ #include "DataCollector.h" #include "FactorHiGHSSettings.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { @@ -93,7 +93,7 @@ void DataCollector::setNorms(double norm1, double maxdiag) { // Printing // -void DataCollector::printTimes(const Log& log) const { +void DataCollector::printTimes(const Logger& logger) const { #if HIPO_TIMING_LEVEL >= 1 std::stringstream log_stream; @@ -289,11 +289,11 @@ void DataCollector::printTimes(const Log& log) const { log_stream << "----------------------------------------------------\n"; #endif - log.print(log_stream); + logger.print(log_stream.str().c_str()); #endif } -void DataCollector::printIter(const Log& log) const { +void DataCollector::printIter(const Logger& logger) const { #ifdef HIPO_COLLECT_EXPENSIVE_DATA std::stringstream log_stream; @@ -319,7 +319,7 @@ void DataCollector::printIter(const Log& log) const { log_stream << sci(iter.M_norm1, 9, 1) << " "; log_stream << sci(iter.M_maxdiag, 9, 1) << "|\n"; } - log.print(log_stream); + logger.print(log_stream.str().c_str()); #endif } diff --git a/highs/ipm/hipo/factorhighs/DataCollector.h b/highs/ipm/hipo/factorhighs/DataCollector.h index 2f67d5752f..e5967787d8 100644 --- a/highs/ipm/hipo/factorhighs/DataCollector.h +++ b/highs/ipm/hipo/factorhighs/DataCollector.h @@ -9,7 +9,7 @@ #include "FactorHiGHSSettings.h" #include "Timing.h" #include "ipm/hipo/auxiliary/IntConfig.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { @@ -77,8 +77,8 @@ class DataCollector { double maxoffD); void setNorms(double norm1, double maxdiag); - void printTimes(const Log& log) const; - void printIter(const Log& log) const; + void printTimes(const Logger& logger) const; + void printIter(const Logger& logger) const; }; } // namespace hipo diff --git a/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp b/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp index 259445cfd0..666703f3cf 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp @@ -5,7 +5,6 @@ #include "ReturnValues.h" #include "Swaps.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" namespace hipo { diff --git a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp index d15c950dd0..c30c44e2b1 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp @@ -5,7 +5,6 @@ #include "ReturnValues.h" #include "Swaps.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" #include "util/HighsRandom.h" namespace hipo { diff --git a/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp b/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp index 8c69cfd97d..3773532f8e 100644 --- a/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp +++ b/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp @@ -3,27 +3,27 @@ #include "Analyse.h" #include "DataCollector.h" #include "Factorise.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { -FHsolver::FHsolver(const Log* log, Int block_size) - : log_{log}, nb_{block_size > 0 ? block_size : default_nb_} { +FHsolver::FHsolver(const Logger* logger, Int block_size) + : logger_{logger}, nb_{block_size > 0 ? block_size : default_nb_} { #ifdef HIPO_COLLECT_EXPENSIVE_DATA - if (log_) - log_->printw( + if (logger_) + logger_->printw( "Running in debug mode: COLLECTING EXPENSIVE FACTORISATION DATA\n"); #endif #if HIPO_TIMING_LEVEL > 0 - if (log_) - log_->printw("Running in debug mode: COLLECTING EXPENSIVE TIMING DATA\n"); + if (logger_) + logger_->printw("Running in debug mode: COLLECTING EXPENSIVE TIMING DATA\n"); #endif } FHsolver::~FHsolver() { - if (log_) { - data_.printTimes(*log_); - data_.printIter(*log_); + if (logger_) { + data_.printTimes(*logger_); + data_.printIter(*logger_); } } @@ -38,14 +38,14 @@ Int FHsolver::analyse(Symbolic& S, const std::vector& rows, const std::vector& ptr, const std::vector& signs, const std::vector& perm) { - Analyse an_obj(rows, ptr, signs, nb_, log_, data_, perm); + Analyse an_obj(rows, ptr, signs, nb_, logger_, data_, perm); return an_obj.run(S); } Int FHsolver::factorise(const Symbolic& S, const std::vector& rows, const std::vector& ptr, const std::vector& vals) { - Factorise fact_obj(S, rows, ptr, vals, regul_, log_, data_, sn_columns_, + Factorise fact_obj(S, rows, ptr, vals, regul_, logger_, data_, sn_columns_, &serial_stack_); return fact_obj.run(N_); } diff --git a/highs/ipm/hipo/factorhighs/FactorHiGHS.h b/highs/ipm/hipo/factorhighs/FactorHiGHS.h index e61141b6dd..79a325db07 100644 --- a/highs/ipm/hipo/factorhighs/FactorHiGHS.h +++ b/highs/ipm/hipo/factorhighs/FactorHiGHS.h @@ -5,7 +5,7 @@ #include "DataCollector.h" #include "Numeric.h" #include "Symbolic.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" /* @@ -37,13 +37,11 @@ Then, the factorization is performed as follows. FH.factorise(S, rows, ptr, val); FH.solve(x); -Printing to screen is achieved using the interface in auxiliary/Log.h. Pass an -object of type Log for normal printing: +Printing to screen is achieved using the interface in auxiliary/Logger.h. ... - Log log; - FHsolver FH(&log); + Logger logger; + FHsolver FH(&logger); ... -Pass an object of type LogHighs for Highs logging. Pass nothing to suppress all logging. To add static regularisation when the pivots are selected, use @@ -58,7 +56,7 @@ input to the constructor. namespace hipo { class FHsolver { - const Log* log_; + const Logger* logger_; DataCollector data_; Regul regul_; Numeric N_; @@ -74,7 +72,7 @@ class FHsolver { public: // Create object and initialise DataCollector - FHsolver(const Log* log = nullptr, Int block_size = default_nb_); + FHsolver(const Logger* logger = nullptr, Int block_size = default_nb_); // Print collected data (if any) and terminate DataCollector ~FHsolver(); diff --git a/highs/ipm/hipo/factorhighs/Factorise.cpp b/highs/ipm/hipo/factorhighs/Factorise.cpp index d03e203d84..0ffeb5ec71 100644 --- a/highs/ipm/hipo/factorhighs/Factorise.cpp +++ b/highs/ipm/hipo/factorhighs/Factorise.cpp @@ -9,7 +9,7 @@ #include "HybridHybridFormatHandler.h" #include "ReturnValues.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "parallel/HighsParallel.h" namespace hipo { @@ -17,13 +17,13 @@ namespace hipo { Factorise::Factorise(const Symbolic& S, const std::vector& rowsM, const std::vector& ptrM, const std::vector& valM, const Regul& regul, - const Log* log, DataCollector& data, + const Logger* logger, DataCollector& data, std::vector>& sn_columns, CliqueStack* stack) : S_{S}, sn_columns_{sn_columns}, regul_{regul}, - log_{log}, + logger_{logger}, data_{data}, stack_{stack} { // Input the symmetric matrix to be factorised in CSC format and the symbolic @@ -33,8 +33,8 @@ Factorise::Factorise(const Symbolic& S, const std::vector& rowsM, n_ = ptrM.size() - 1; if (n_ != S_.size()) { - if (log_) - log_->printDevInfo( + if (logger_) + logger_->printInfo( "Matrix provided to Factorise has size incompatible with symbolic " "object.\n"); return; @@ -134,8 +134,8 @@ void Factorise::processSupernode(Int sn) { if (serial) { bool reallocation = false; clique_ptr = stack_->setup(S_.cliqueSize(sn), reallocation); - if (reallocation && log_) - log_->printDevInfo("Reallocation of CliqueStack\n"); + if (reallocation && logger_) + logger_->printInfo("Reallocation of CliqueStack\n"); } // initialise the format handler @@ -185,7 +185,7 @@ void Factorise::processSupernode(Int sn) { child_clique = schur_contribution_[child_sn].data(); if (!child_clique) { - if (log_) log_->printDevInfo("Missing child supernode contribution\n"); + if (logger_) logger_->printInfo("Missing child supernode contribution\n"); flag_stop_.store(true, std::memory_order_relaxed); return; } @@ -263,10 +263,10 @@ void Factorise::processSupernode(Int sn) { if (Int flag = FH->denseFactorise(reg_thresh)) { flag_stop_.store(true, std::memory_order_relaxed); - if (log_ && flag == kRetInvalidInput) - log_->printDevInfo("DenseFact: invalid input\n"); - else if (log_ && flag == kRetInvalidPivot) - log_->printDevInfo("DenseFact: invalid pivot\n"); + if (logger_ && flag == kRetInvalidInput) + logger_->printInfo("DenseFact: invalid input\n"); + else if (logger_ && flag == kRetInvalidPivot) + logger_->printInfo("DenseFact: invalid pivot\n"); return; } diff --git a/highs/ipm/hipo/factorhighs/Factorise.h b/highs/ipm/hipo/factorhighs/Factorise.h index d7887223c9..0a673b5c50 100644 --- a/highs/ipm/hipo/factorhighs/Factorise.h +++ b/highs/ipm/hipo/factorhighs/Factorise.h @@ -7,7 +7,7 @@ #include "Numeric.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { @@ -62,7 +62,7 @@ class Factorise { // flag to stop computation std::atomic flag_stop_{false}; - const Log* log_; + const Logger* logger_; DataCollector& data_; CliqueStack* stack_; @@ -74,7 +74,7 @@ class Factorise { public: Factorise(const Symbolic& S, const std::vector& rowsM, const std::vector& ptrM, const std::vector& valM, - const Regul& regul, const Log* log, DataCollector& data, + const Regul& regul, const Logger* logger, DataCollector& data, std::vector>& sn_columns, CliqueStack* stack); bool run(Numeric& num); diff --git a/highs/ipm/hipo/factorhighs/Numeric.cpp b/highs/ipm/hipo/factorhighs/Numeric.cpp index 963a5fd933..5872d5d9d3 100644 --- a/highs/ipm/hipo/factorhighs/Numeric.cpp +++ b/highs/ipm/hipo/factorhighs/Numeric.cpp @@ -6,7 +6,6 @@ #include "ReturnValues.h" #include "Timing.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" #include "ipm/hipo/auxiliary/VectorOperations.h" #include "util/HighsCDouble.h" #include "util/HighsRandom.h" diff --git a/highs/ipm/hipo/factorhighs/Symbolic.cpp b/highs/ipm/hipo/factorhighs/Symbolic.cpp index b4626714af..163154281d 100644 --- a/highs/ipm/hipo/factorhighs/Symbolic.cpp +++ b/highs/ipm/hipo/factorhighs/Symbolic.cpp @@ -3,7 +3,6 @@ #include #include "FactorHiGHSSettings.h" -#include "ipm/hipo/auxiliary/Log.h" namespace hipo { @@ -64,7 +63,7 @@ static std::string memoryString(double mem) { return ss.str(); } -void Symbolic::print(const Log& log, bool verbose) const { +void Symbolic::print(const Logger& logger, bool verbose) const { std::stringstream log_stream; log_stream << "\nFactorisation statistics\n"; log_stream << textline("Size:") << sci(n_, 0, 2) << '\n'; @@ -95,7 +94,7 @@ void Symbolic::print(const Log& log, bool verbose) const { log_stream << textline("Sn avg size:") << sci((double)n_ / sn_, 0, 1) << '\n'; } - log.print(log_stream); + logger.print(log_stream.str().c_str()); } } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/Symbolic.h b/highs/ipm/hipo/factorhighs/Symbolic.h index 2fd9e1faad..1c1dff63d6 100644 --- a/highs/ipm/hipo/factorhighs/Symbolic.h +++ b/highs/ipm/hipo/factorhighs/Symbolic.h @@ -4,7 +4,7 @@ #include #include "ipm/hipo/auxiliary/IntConfig.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" namespace hipo { @@ -133,7 +133,7 @@ class Symbolic { const std::vector& snStart() const; const std::vector& pivotSign() const; - void print(const Log& log, bool verbose = false) const; + void print(const Logger& logger, bool verbose = false) const; }; // Explanation of relative indices: diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 61d946a722..561b5d1c7c 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -5,7 +5,7 @@ #include "Status.h" #include "amd/amd.h" #include "ipm/hipo/auxiliary/Auxiliary.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "metis/metis.h" #include "parallel/HighsParallel.h" #include "rcm/rcm.h" @@ -15,14 +15,14 @@ namespace hipo { FactorHiGHSSolver::FactorHiGHSSolver(KktMatrix& kkt, Options& options, const Model& model, const Regularisation& regul, Info& info, - IpmData& record, const LogHighs& log) - : FH_(&log, options.block_size), + IpmData& record, const Logger& logger) + : FH_(&logger, options.block_size), S_{}, kkt_{kkt}, regul_{regul}, info_{info}, data_{record}, - log_{log}, + logger_{logger}, model_{model}, options_{options} {} @@ -48,7 +48,7 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) { std::vector pivot_signs(n + m, -1); for (Int i = 0; i < m; ++i) pivot_signs[n + i] = 1; - log_.printDevInfo("Performing AS analyse phase\n"); + logger_.printInfo("Performing AS analyse phase\n"); Clock clock; Int status = chooseOrdering(kkt_.rowsAS, kkt_.ptrAS, pivot_signs, S, @@ -68,7 +68,7 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S) { // create vector of signs of pivots std::vector pivot_signs(model_.A().num_row_, 1); - log_.printDevInfo("Performing NE analyse phase\n"); + logger_.printInfo("Performing NE analyse phase\n"); Clock clock; Int status = chooseOrdering(kkt_.rowsNE, kkt_.ptrNE, pivot_signs, S, @@ -184,16 +184,16 @@ Int FactorHiGHSSolver::setup() { std::stringstream log_stream; log_stream << textline("Analyse time:") << fix(clock.stop(), 0, 2) << '\n'; - log_.print(log_stream); + logger_.print(log_stream.str().c_str()); - S_.print(log_, log_.debug(1)); + S_.print(logger_, logger_.debug(1)); // Warn about large memory consumption if (S_.storage() > kLargeStorageGB * 1024 * 1024 * 1024) { - log_.printw("Large amount of memory required\n"); + logger_.printw("Large amount of memory required\n"); } - log_.print("\n"); + logger_.print("\n"); return kStatusOk; } @@ -214,13 +214,13 @@ Int FactorHiGHSSolver::chooseNla() { if (expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { failure_NE = true; - log_.printDevInfo("NE skipped\n"); + logger_.printInfo("NE skipped\n"); } else { Int status = kkt_.buildNEstructure(); if (status) { failure_NE = true; if (status == kStatusOverflow) { - log_.printDevInfo("Integer overflow forming NE matrix\n"); + logger_.printInfo("Integer overflow forming NE matrix\n"); overflow_NE = true; } return; @@ -242,13 +242,13 @@ Int FactorHiGHSSolver::chooseNla() { if (expect_NE_much_cheaper && can_skip_AS) { failure_AS = true; - log_.printDevInfo("AS skipped\n"); + logger_.printInfo("AS skipped\n"); } else { Int AS_status = kkt_.buildASstructure(); if (!AS_status) AS_status = analyseAS(symb_AS); if (AS_status) failure_AS = true; if (AS_status == kStatusOverflow) { - log_.printDevInfo("Integer overflow forming AS matrix\n"); + logger_.printInfo("Integer overflow forming AS matrix\n"); overflow_AS = true; } @@ -291,7 +291,7 @@ Int FactorHiGHSSolver::chooseNla() { else status = kStatusErrorAnalyse; - log_.printe("Both NE and AS failed analyse phase\n"); + logger_.printe("Both NE and AS failed analyse phase\n"); } else { // Total number of operations, given by dense flops and sparse indexing // operations, weighted with an empirical factor @@ -319,7 +319,7 @@ Int FactorHiGHSSolver::chooseNla() { } } - log_.print(log_stream); + logger_.print(log_stream.str().c_str()); if (status == kStatusOk) { if (options_.nla == kHipoAugmentedString) { @@ -358,7 +358,7 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, if (nla == "NE") { if (ptr.back() >= kkt_.NE_nz_limit.load(std::memory_order_relaxed)) { - log_.printDevInfo("NE interrupted\n"); + logger_.printInfo("NE interrupted\n"); return kStatusErrorAnalyse; } } @@ -383,22 +383,22 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // set logging of Metis depending on debug level options[METIS_OPTION_DBGLVL] = 0; - if (log_.debug(2)) + if (logger_.debug(2)) options[METIS_OPTION_DBGLVL] = METIS_DBG_INFO | METIS_DBG_COARSEN; // no2hop improves the quality of ordering in general options[METIS_OPTION_NO2HOP] = 1; - log_.printDevInfo("Running Metis\n"); + logger_.printInfo("Running Metis\n"); std::vector iperm(n); Int status = Highs_METIS_NodeND(&n, full_ptr.data(), full_rows.data(), NULL, options, permutations[i].data(), iperm.data()); - log_.printDevInfo("Metis done\n"); + logger_.printInfo("Metis done\n"); if (status != METIS_OK) { - log_.printDevInfo("Error with Metis\n"); + logger_.printInfo("Error with Metis\n"); failure[i] = true; } } else if (orderings_to_try[i] == kHipoAmdString) { @@ -406,23 +406,23 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, Highs_amd_defaults(control); double info[AMD_INFO]; - log_.printDevInfo("Running AMD\n"); + logger_.printInfo("Running AMD\n"); Int status = Highs_amd_order(n, full_ptr.data(), full_rows.data(), permutations[i].data(), control, info); - log_.printDevInfo("AMD done\n"); + logger_.printInfo("AMD done\n"); if (status != AMD_OK) { - log_.printDevInfo("Error with AMD\n"); + logger_.printInfo("Error with AMD\n"); failure[i] = true; } } else if (orderings_to_try[i] == kHipoRcmString) { - log_.printDevInfo("Running RCM\n"); + logger_.printInfo("Running RCM\n"); Int status = Highs_genrcm(n, full_ptr.back(), full_ptr.data(), full_rows.data(), permutations[i].data()); - log_.printDevInfo("RCM done\n"); + logger_.printInfo("RCM done\n"); if (status != 0) { - log_.printDevInfo("Error with RCM\n"); + logger_.printInfo("Error with RCM\n"); failure[i] = true; } } else { @@ -433,9 +433,9 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, failure[i] = FH_.analyse(symbolics[i], rows, ptr, signs, permutations[i]); - if (failure[i] && log_.debug(2)) { - log_.print("Failed symbolic:"); - symbolics[i].print(log_, true); + if (failure[i] && logger_.debug(2)) { + logger_.print("Failed symbolic:"); + symbolics[i].print(logger_, true); } }; @@ -518,7 +518,7 @@ Int FactorHiGHSSolver::setNla() { std::stringstream log_stream; if (options_.nla == kHipoNormalEqString && model_.nonSeparableQp()) { - log_.printw("Normal equations not available for non-separable QP\n"); + logger_.printw("Normal equations not available for non-separable QP\n"); options_.nla = kHighsChooseString; } @@ -526,10 +526,10 @@ Int FactorHiGHSSolver::setNla() { Int status = kkt_.buildASstructure(); if (!status) status = analyseAS(S_); if (status == kStatusOverflow) { - log_.printe("AS requested, integer overflow\n"); + logger_.printe("AS requested, integer overflow\n"); return kStatusOverflow; } else if (status) { - log_.printe("AS requested, failed analyse phase\n"); + logger_.printe("AS requested, failed analyse phase\n"); return kStatusErrorAnalyse; } log_stream << textline("Newton system:") << "AS requested\n"; @@ -538,10 +538,10 @@ Int FactorHiGHSSolver::setNla() { Int status = kkt_.buildNEstructure(); if (!status) status = analyseNE(S_); if (status == kStatusOverflow) { - log_.printe("NE requested, integer overflow\n"); + logger_.printe("NE requested, integer overflow\n"); return kStatusOverflow; } else if (status) { - log_.printe("NE requested, failed analyse phase\n"); + logger_.printe("NE requested, failed analyse phase\n"); return kStatusErrorAnalyse; } log_stream << textline("Newton system:") << "NE requested\n"; @@ -552,12 +552,12 @@ Int FactorHiGHSSolver::setNla() { } else assert(1 == 0); - if (log_.debug(1)) + if (logger_.debug(1)) log_stream << textline("Ordering:") << (options_.nla == kHipoAugmentedString ? ordering_AS_ : ordering_NE_) << '\n'; - log_.print(log_stream); + logger_.print(log_stream.str().c_str()); kkt_.iperm = S_.iperm(); @@ -641,7 +641,7 @@ void FactorHiGHSSolver::setParallel() { } else assert(1 == 0); - log_.print(log_stream); + logger_.print(log_stream.str().c_str()); S_.setParallel(parallel_tree, parallel_node); } diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h index 779758b363..651147aadb 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h @@ -26,7 +26,7 @@ class FactorHiGHSSolver : public LinearSolver { const Regularisation& regul_; Info& info_; IpmData& data_; - const LogHighs& log_; + const Logger& logger_; const Model& model_; Options& options_; @@ -45,7 +45,7 @@ class FactorHiGHSSolver : public LinearSolver { public: FactorHiGHSSolver(KktMatrix& kkt, Options& options, const Model& model, const Regularisation& regul, Info& info, IpmData& record, - const LogHighs& log); + const Logger& logger); // Override functions Int factorAS(const std::vector& scaling) override; diff --git a/highs/ipm/hipo/ipm/KktMatrix.cpp b/highs/ipm/hipo/ipm/KktMatrix.cpp index 1926c4c337..304ea874ff 100644 --- a/highs/ipm/hipo/ipm/KktMatrix.cpp +++ b/highs/ipm/hipo/ipm/KktMatrix.cpp @@ -5,8 +5,8 @@ namespace hipo { KktMatrix::KktMatrix(const Model& m, const Regularisation& r, Info& i, - const LogHighs& l) - : model{m}, regul{r}, info{i}, log{l} {} + const Logger& l) + : model{m}, regul{r}, info{i}, logger{l} {} Int KktMatrix::buildASstructure() { // Build lower triangular structure of the augmented system. @@ -20,7 +20,7 @@ Int KktMatrix::buildASstructure() { const Int nzA = A.numNz(); const Int nzQ = Q.numNz(); - log.printDevInfo("Building AS structure\n"); + logger.printInfo("Building AS structure\n"); const Int nzBlock11 = model.qp() ? nzQ : n; @@ -101,7 +101,7 @@ Int KktMatrix::buildNEstructure() { const Int m = A.num_row_; const Int nzA = A.numNz(); - log.printDevInfo("Building NE structure\n"); + logger.printInfo("Building NE structure\n"); // create partial row-wise representation without values, and array or // corresponding indices between cw and rw representation diff --git a/highs/ipm/hipo/ipm/KktMatrix.h b/highs/ipm/hipo/ipm/KktMatrix.h index 6d64c89dbd..c83ea75402 100644 --- a/highs/ipm/hipo/ipm/KktMatrix.h +++ b/highs/ipm/hipo/ipm/KktMatrix.h @@ -27,10 +27,10 @@ struct KktMatrix { const Model& model; const Regularisation& regul; Info& info; - const LogHighs& log; + const Logger& logger; KktMatrix(const Model& model, const Regularisation& regul, Info& info, - const LogHighs& log); + const Logger& logger); Int buildASstructure(); Int buildASvalues(const std::vector& scaling); diff --git a/highs/ipm/hipo/ipm/LogHighs.cpp b/highs/ipm/hipo/ipm/LogHighs.cpp deleted file mode 100644 index 98417ed19b..0000000000 --- a/highs/ipm/hipo/ipm/LogHighs.cpp +++ /dev/null @@ -1,71 +0,0 @@ -#include "LogHighs.h" - -namespace hipo { - -void LogHighs::setOptions(const HighsLogOptions* log_options) { - log_options_ = log_options; -} - -bool LogHighs::debug(Int level) const { - // Return true if level agrees with log_dev_level - - if (log_options_ && log_options_->log_dev_level) { - if (*log_options_->log_dev_level == kHighsLogDevLevelInfo) - return level == 1; - - if (*log_options_->log_dev_level == kHighsLogDevLevelDetailed) - return level == 1 || level == 2; - - if (*log_options_->log_dev_level == kHighsLogDevLevelVerbose) - return level == 1 || level == 2 || level == 3; - } - return false; -} - -void LogHighs::print(std::stringstream& ss) const { - if (log_options_) - highsLogUser(*log_options_, HighsLogType::kInfo, "%s", ss.str().c_str()); -} -void LogHighs::printw(std::stringstream& ss) const { - if (log_options_) - highsLogUser(*log_options_, HighsLogType::kWarning, "%s", ss.str().c_str()); -} -void LogHighs::printe(std::stringstream& ss) const { - if (log_options_) - highsLogUser(*log_options_, HighsLogType::kError, "%s", ss.str().c_str()); -} -void LogHighs::print(const char* c) const { - if (log_options_) highsLogUser(*log_options_, HighsLogType::kInfo, "%s", c); -} -void LogHighs::printw(const char* c) const { - if (log_options_) - highsLogUser(*log_options_, HighsLogType::kWarning, "%s", c); -} -void LogHighs::printe(const char* c) const { - if (log_options_) highsLogUser(*log_options_, HighsLogType::kError, "%s", c); -} - -void LogHighs::printDevInfo(std::stringstream& ss) const { - if (log_options_) - highsLogDev(*log_options_, HighsLogType::kInfo, "%s", ss.str().c_str()); -} -void LogHighs::printDevDetailed(std::stringstream& ss) const { - if (log_options_) - highsLogDev(*log_options_, HighsLogType::kDetailed, "%s", ss.str().c_str()); -} -void LogHighs::printDevVerbose(std::stringstream& ss) const { - if (log_options_) - highsLogDev(*log_options_, HighsLogType::kVerbose, "%s", ss.str().c_str()); -} -void LogHighs::printDevInfo(const char* c) const { - if (log_options_) highsLogDev(*log_options_, HighsLogType::kInfo, "%s", c); -} -void LogHighs::printDevDetailed(const char* c) const { - if (log_options_) - highsLogDev(*log_options_, HighsLogType::kDetailed, "%s", c); -} -void LogHighs::printDevVerbose(const char* c) const { - if (log_options_) highsLogDev(*log_options_, HighsLogType::kVerbose, "%s", c); -} - -} // namespace hipo diff --git a/highs/ipm/hipo/ipm/LogHighs.h b/highs/ipm/hipo/ipm/LogHighs.h deleted file mode 100644 index eedb6a73e5..0000000000 --- a/highs/ipm/hipo/ipm/LogHighs.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef HIPO_LOG_HIGHS_H -#define HIPO_LOG_HIGHS_H - -#include "ipm/hipo/auxiliary/Log.h" - -// Class for Highs logging -// Use setOptions to pass the log_options. - -namespace hipo { - -class LogHighs : public Log { - const HighsLogOptions* log_options_; - - public: - void setOptions(const HighsLogOptions* log_options); - bool debug(Int level) const override; - - void print(std::stringstream& ss) const override; - void printw(std::stringstream& ss) const override; - void printe(std::stringstream& ss) const override; - void print(const char* c) const override; - void printw(const char* c) const override; - void printe(const char* c) const override; - - void printDevInfo(std::stringstream& ss) const override; - void printDevDetailed(std::stringstream& ss) const override; - void printDevVerbose(std::stringstream& ss) const override; - void printDevInfo(const char* c) const override; - void printDevDetailed(const char* c) const override; - void printDevVerbose(const char* c) const override; -}; -} // namespace hipo -#endif \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index e61acf3655..4cb73407a8 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -3,7 +3,7 @@ #include "Parameters.h" #include "Status.h" #include "ipm/IpxWrapper.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "model/HighsHessianUtils.h" namespace hipo { @@ -177,7 +177,7 @@ void Model::computeNorms() { } } -void Model::print(const LogHighs& log) const { +void Model::print(const Logger& logger) const { std::stringstream log_stream; log_stream << textline("Rows:") << sci(m_, 0, 1) << '\n'; @@ -308,7 +308,7 @@ void Model::print(const LogHighs& log) const { log_stream << textline("Range of A:") << "[" << sci(Amin, 5, 1) << ", " << sci(Amax, 5, 1) << "]\n"; - if (log.debug(1)) { + if (logger.debug(1)) { log_stream << textline("Inf-norm rows:") << "[" << sci(norm_inf_row_min, 5, 1) << ", " << sci(norm_inf_row_max, 5, 1) << "]\n"; @@ -343,7 +343,7 @@ void Model::print(const LogHighs& log) const { log_stream << textline("Scaling coefficients:") << "[" << sci(scalemin, 5, 1) << ", " << sci(scalemax, 5, 1) << "]\n"; - if (log.debug(1)) { + if (logger.debug(1)) { preprocessor_.print(log_stream); log_stream << textline("Free variables:") << integer(free_vars) << '\n'; @@ -353,7 +353,7 @@ void Model::print(const LogHighs& log) const { << "]\n"; } - log.print(log_stream); + logger.print(log_stream.str().c_str()); } Int Model::loadIntoIpx(ipx::LpSolver& lps) const { diff --git a/highs/ipm/hipo/ipm/Model.h b/highs/ipm/hipo/ipm/Model.h index e36541e50c..1d3c01120e 100644 --- a/highs/ipm/hipo/ipm/Model.h +++ b/highs/ipm/hipo/ipm/Model.h @@ -5,7 +5,7 @@ #include #include -#include "LogHighs.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "PreProcess.h" #include "ipm/hipo/auxiliary/IntConfig.h" #include "ipm/ipx/lp_solver.h" @@ -74,7 +74,7 @@ class Model { Int init(const HighsLp& lp, const HighsHessian& Q); // Print information of model - void print(const LogHighs& log) const; + void print(const Logger& logger) const; void printDense() const; diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index 4334825ba6..1f2400c62c 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -5,7 +5,7 @@ #include #include "ipm/IpxWrapper.h" -#include "ipm/hipo/auxiliary/Log.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "lp_data/HighsSolution.h" #include "parallel/HighsParallel.h" @@ -13,7 +13,7 @@ namespace hipo { Int Solver::load(const HighsLp& lp, const HighsHessian& Q) { if (model_.init(lp, Q)) { - logH_.printDevInfo("Error with model\n"); + logger_.printInfo("Error with model\n"); return kStatusBadModel; } return kStatusOk; @@ -27,7 +27,7 @@ void Solver::setOptions(const HighsOptions& highs_options) { options_.log_options = &highs_options.log_options; // Debug option is already considered through log_options.log_dev_level in - // hipo::LogHighs::debug + // hipo::Logger::debug options_.timeless_log = highs_options.timeless_log; options_.feasibility_tol = @@ -60,7 +60,7 @@ void Solver::setOptions(const HighsOptions& highs_options) { void Solver::resetOptions() { options_ = options_orig_; - if (options_.display) logH_.setOptions(options_.log_options); + if (options_.display) logger_.setOptions(options_.log_options); control_.setOptions(options_); } void Solver::setCallback(HighsCallback& callback) { @@ -136,7 +136,7 @@ bool Solver::initialise() { start_time_ = control_.elapsed(); - kkt_.reset(new KktMatrix(model_, regul_, info_, logH_)); + kkt_.reset(new KktMatrix(model_, regul_, info_, logger_)); if (!kkt_) { info_.status = kStatusError; return true; @@ -144,7 +144,7 @@ bool Solver::initialise() { // initialise linear solver LS_.reset(new FactorHiGHSSolver(*kkt_, options_, model_, regul_, info_, - it_->data, logH_)); + it_->data, logger_)); if (Int status = LS_->setup()) { info_.status = (Status)status; return true; @@ -259,7 +259,7 @@ bool Solver::prepareIpx() { Int load_status = model_.loadIntoIpx(ipx_lps_); if (load_status) { - logH_.printDevInfo("Error loading model into IPX\n"); + logger_.printInfo("Error loading model into IPX\n"); return true; } @@ -271,7 +271,7 @@ bool Solver::prepareIpx() { zu.data()); if (start_point_status) { - logH_.printDevInfo("Error loading starting point into IPX\n"); + logger_.printInfo("Error loading starting point into IPX\n"); return true; } @@ -282,9 +282,9 @@ void Solver::refineWithIpx() { if (checkInterrupt()) return; if (statusNeedsRefinement() && refinementIsOn()) { - logH_.print("\nRestarting with IPX\n"); + logger_.print("\nRestarting with IPX\n"); } else if (statusAllowsCrossover() && crossoverIsOn()) { - logH_.print("\nRunning crossover with IPX\n"); + logger_.print("\nRunning crossover with IPX\n"); } else { return; } @@ -306,7 +306,7 @@ void Solver::refineWithIpx() { log_stream << ", crossover " << ipx::StatusString(info_.ipx_info.status_crossover); log_stream << '\n'; - logH_.print(log_stream); + logger_.print(log_stream.str().c_str()); if (info_.ipx_info.status_crossover == IPX_STATUS_optimal) { info_.status = kStatusBasic; @@ -319,11 +319,11 @@ bool Solver::solveNewtonSystem(NewtonDir& delta) { // Check for NaN of Inf if (it_->isDirNan(delta)) { - logH_.printDevInfo("Direction is nan\n"); + logger_.printInfo("Direction is nan\n"); info_.status = kStatusError; return true; } else if (it_->isDirInf(delta)) { - logH_.printDevInfo("Direction is inf\n"); + logger_.printInfo("Direction is inf\n"); info_.status = kStatusError; return true; } @@ -343,7 +343,7 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) { // factorise normal equations, if not yet done if (!LS_->valid_) { if (Int status = LS_->factorNE(theta_inv)) { - logH_.printe("Error while factorising normal equations\n"); + logger_.printe("Error while factorising normal equations\n"); info_.status = (Status)status; return true; } @@ -352,7 +352,7 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) { // solve with normal equations if (Int status = LS_->solveNE(res8, delta.y)) { - logH_.printe("Error while solving normal equations\n"); + logger_.printe("Error while solving normal equations\n"); info_.status = (Status)status; return true; } @@ -377,7 +377,7 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) { // factorise augmented system, if not yet done if (!LS_->valid_) { if (Int status = LS_->factorAS(theta_inv)) { - logH_.printe("Error while factorising augmented system\n"); + logger_.printe("Error while factorising augmented system\n"); info_.status = (Status)status; return true; } @@ -386,7 +386,7 @@ bool Solver::solve2x2(NewtonDir& delta, const Residuals& rhs) { // solve with augmented system if (Int status = LS_->solveAS(res7, rhs.r1, delta.x, delta.y)) { - logH_.printe("Error while solving augmented system\n"); + logger_.printe("Error while solving augmented system\n"); info_.status = (Status)status; return true; } @@ -613,13 +613,13 @@ bool Solver::startingPoint() { // factorise A*A^T if (Int status = LS_->factorNE(empty_scaling)) { - logH_.printe("Error while factorising normal equations\n"); + logger_.printe("Error while factorising normal equations\n"); info_.status = (Status)status; return true; } if (Int status = LS_->solveNE(y, temp_m)) { - logH_.printe("Error while solving normal equations\n"); + logger_.printe("Error while solving normal equations\n"); info_.status = (Status)status; return true; } @@ -630,7 +630,7 @@ bool Solver::startingPoint() { // [ A 0 ] [ dx] = [ b ] if (Int status = LS_->factorAS(empty_scaling)) { - logH_.printe("Error while factorising augmented system\n"); + logger_.printe("Error while factorising augmented system\n"); info_.status = (Status)status; return true; } @@ -639,7 +639,7 @@ bool Solver::startingPoint() { for (Int i = 0; i < n_; ++i) rhs_x[i] = -x[i]; std::vector lhs_x(n_); if (Int status = LS_->solveAS(rhs_x, y, lhs_x, temp_m)) { - logH_.printe("Error while solving augmented system\n"); + logger_.printe("Error while solving augmented system\n"); info_.status = (Status)status; return true; } @@ -698,7 +698,7 @@ bool Solver::startingPoint() { model_.A().alphaProductPlusY(1.0, cPlusQx, temp_m); if (Int status = LS_->solveNE(temp_m, y)) { - logH_.printe("Error while solvingF normal equations\n"); + logger_.printe("Error while solvingF normal equations\n"); info_.status = (Status)status; return true; } @@ -712,7 +712,7 @@ bool Solver::startingPoint() { std::vector lhs_x(n_); if (Int status = LS_->solveAS(cPlusQx, rhs_y, lhs_x, y)) { - logH_.printe("Error while solving augmented system\n"); + logger_.printe("Error while solving augmented system\n"); info_.status = (Status)status; return true; } @@ -992,11 +992,11 @@ void Solver::bestWeight(const NewtonDir& delta, const NewtonDir& corrector, bool Solver::checkIterate() { // Check that iterate is not NaN or Inf if (it_->isNan()) { - logH_.printDevInfo("\nIterate is nan\n"); + logger_.printInfo("\nIterate is nan\n"); info_.status = kStatusError; return true; } else if (it_->isInf()) { - logH_.printDevInfo("\nIterate is inf\n"); + logger_.printInfo("\nIterate is inf\n"); info_.status = kStatusError; return true; } @@ -1007,7 +1007,7 @@ bool Solver::checkIterate() { (model_.hasLb(i) && it_->zl[i] < 0) || (model_.hasUb(i) && it_->xu[i] < 0) || (model_.hasUb(i) && it_->zu[i] < 0)) { - logH_.printDevInfo("\nIterate has negative component\n"); + logger_.printInfo("\nIterate has negative component\n"); return true; } } @@ -1018,7 +1018,7 @@ bool Solver::checkIterate() { bool Solver::checkStagnation() { std::stringstream log_stream; bool stagnation = it_->stagnation(log_stream); - logH_.printDevInfo(log_stream); + logger_.printInfo(log_stream.str().c_str()); return stagnation; } @@ -1046,20 +1046,20 @@ bool Solver::checkBadIter() { if (pobj_is_larger) { // problem is likely to be primal unbounded, i.e. dual infeasible - logH_.print("=== Dual infeasible\n"); + logger_.print("=== Dual infeasible\n"); info_.status = kStatusDualInfeasible; terminate = true; } else if (dobj_is_larger) { // problem is likely to be dual unbounded, i.e. primal infeasible - logH_.print("=== Primal infeasible\n"); + logger_.print("=== Primal infeasible\n"); info_.status = kStatusPrimalInfeasible; terminate = true; } else if (stagnation) { // stagnation detected, solution may still be good for highs kktCheck if (info_.status != kStatusPDFeas && checkTerminationKkt()) { - logH_.printw( + logger_.printw( "HiPO stagnated but HiGHS considers the solution acceptable\n"); - logH_.print("=== Primal-dual feasible point found\n"); + logger_.print("=== Primal-dual feasible point found\n"); info_.status = kStatusPDFeas; } else info_.status = kStatusNoProgress; @@ -1080,12 +1080,12 @@ bool Solver::checkTermination() { if (feasible && optimal) { if (crossoverIsOn()) { if (info_.status != kStatusPDFeas) - logH_.print("=== Primal-dual feasible point found\n"); + logger_.print("=== Primal-dual feasible point found\n"); info_.status = kStatusPDFeas; bool ready_for_crossover = it_->infeasAfterDropping() < options_.crossover_tol; if (ready_for_crossover) { - logH_.print("=== Ready for crossover\n"); + logger_.print("=== Ready for crossover\n"); terminate = true; } @@ -1093,7 +1093,7 @@ bool Solver::checkTermination() { terminate = model_.qp() ? true : checkTerminationKkt(); if (terminate) { assert(info_.status != kStatusPDFeas); - logH_.print("=== Primal-dual feasible point found\n"); + logger_.print("=== Primal-dual feasible point found\n"); info_.status = kStatusPDFeas; } } @@ -1105,14 +1105,14 @@ bool Solver::checkTermination() { bool Solver::checkTerminationKkt() { assert(!model_.qp()); - logH_.printDevInfo("Solution may be optimal, perform kktCheck\n"); + logger_.printInfo("Solution may be optimal, perform kktCheck\n"); HighsModelStatus model_status = HighsModelStatus::kOptimal; HighsInfo highs_info; HighsSolution highs_solution; // Allow kktCheck to print only in debug mode (this is a copy of highs // options, not the original) - Hoptions_.output_flag = logH_.debug(2); + Hoptions_.output_flag = logger_.debug(2); if (!model_.lpOrig()) return false; @@ -1125,10 +1125,10 @@ bool Solver::checkTerminationKkt() { Hoptions_, "During HiPO solve"); if (model_status == HighsModelStatus::kOptimal) { - logH_.printDevInfo("Check successfull\n"); + logger_.printInfo("Check successfull\n"); return true; } else - logH_.printDevInfo("Check failed\n"); + logger_.printInfo("Check failed\n"); return false; } @@ -1145,58 +1145,37 @@ bool Solver::checkInterrupt() { void Solver::printHeader() const { if (iter_ % 20 == 0) { - std::stringstream log_stream; - log_stream << " iter primal obj dual obj" - << " pinf dinf gap"; - - if (!options_.timeless_log) log_stream << " time"; - - if (logH_.debug(1)) { - log_stream << " alpha p/d sigma af/co cor solv" - << " static reg p/d minT maxT" - << " (xj * zj / mu)_range_&_num max_res"; + logger_.print( + " iter primal obj dual obj pinf dinf " + "gap"); + if (!options_.timeless_log) logger_.print(" time"); + if (logger_.debug(1)) { + logger_.print( + " alpha p/d sigma af/co cor solv static reg p/d minT " + " maxT (xj * zj / mu)_range_&_num max_res"); } - - log_stream << "\n"; - logH_.print(log_stream); + logger_.print("\n"); } } void Solver::printOutput() const { printHeader(); - std::stringstream log_stream; - log_stream << integer(iter_, 5); - log_stream << " " << sci(it_->pobj, 16, 8); - log_stream << " " << sci(it_->dobj, 16, 8); - log_stream << " " << sci(it_->pinf, 10, 2); - log_stream << " " << sci(it_->dinf, 10, 2); - log_stream << " " << sci(it_->pdgap, 9, 2); - if (!options_.timeless_log) - log_stream << " " << fix(control_.elapsed(), 7, 1); - - if (logH_.debug(1)) { + logger_.print("%5d %16.8e %16.8e %10.2e %10.2e %9.2e", iter_, it_->pobj, + it_->dobj, it_->pinf, it_->dinf, it_->pdgap); + if (!options_.timeless_log) logger_.print(" %7.1f", control_.elapsed()); + if (logger_.debug(1)) { const IpmIterData& data = it_->data.back(); - - log_stream << " " << fix(alpha_primal_, 6, 2); - log_stream << " " << fix(alpha_dual_, 6, 2); - log_stream << " " << fix(data.sigma_aff, 6, 2); - log_stream << " " << fix(data.sigma, 6, 2); - log_stream << " " << integer(data.correctors, 5); - log_stream << " " << integer(data.num_solves, 5); - log_stream << " " << sci(regul_.primal, 8, 1); - log_stream << " " << sci(regul_.dual, 8, 1); - log_stream << " " << sci(data.min_theta, 8, 1); - log_stream << " " << sci(data.max_theta, 8, 1); - log_stream << " " << sci(data.min_prod, 8, 1); - log_stream << " " << sci(data.max_prod, 8, 1); - log_stream << " " << integer(data.num_small_prod, 4); - log_stream << " " << integer(data.num_large_prod, 4); - log_stream << " " << sci(data.omega, 9, 1); + logger_.print( + " %6.2f %6.2f %6.2f %6.2f %5d %5d %8.1e %8.1e %8.1e %8.1e %8.1e %8.1e " + "%4d " + "%4d %9.1e", + alpha_primal_, alpha_dual_, data.sigma_aff, data.sigma, data.correctors, + data.num_solves, regul_.primal, regul_.dual, data.min_theta, + data.max_theta, data.min_prod, data.max_prod, data.num_small_prod, + data.num_large_prod, data.omega); } - - log_stream << "\n"; - logH_.print(log_stream); + logger_.print("\n"); } void Solver::printInfo() const { @@ -1210,10 +1189,10 @@ void Solver::printInfo() const { log_stream << textline("Threads:") << highs::parallel::num_threads() << '\n'; - logH_.print(log_stream); + logger_.print(log_stream.str().c_str()); // print information about model - model_.print(logH_); + model_.print(logger_); } void Solver::printSummary() const { @@ -1230,7 +1209,7 @@ void Solver::printSummary() const { log_stream << textline("IPX iterations:") << integer(info_.ipx_info.iter) << "\n"; - if (logH_.debug(1)) { + if (logger_.debug(1)) { log_stream << textline("Correctors:") << integer(info_.correctors) << '\n'; log_stream << textline("Factorisations:") << integer(info_.factor_number) << '\n'; @@ -1262,7 +1241,7 @@ void Solver::printSummary() const { } } - logH_.print(log_stream); + logger_.print(log_stream.str().c_str()); } void Solver::getOriginalDims(Int& num_row, Int& num_col) const { diff --git a/highs/ipm/hipo/ipm/Solver.h b/highs/ipm/hipo/ipm/Solver.h index a59789610f..0848667f0d 100644 --- a/highs/ipm/hipo/ipm/Solver.h +++ b/highs/ipm/hipo/ipm/Solver.h @@ -8,7 +8,7 @@ #include "Info.h" #include "Iterate.h" #include "LinearSolver.h" -#include "LogHighs.h" +#include "ipm/hipo/auxiliary/Logger.h" #include "Model.h" #include "Options.h" #include "Parameters.h" @@ -66,7 +66,7 @@ class Solver { ipx::LpSolver ipx_lps_; // Interface to Highs logging - LogHighs logH_; + Logger logger_; double start_time_; From ed5d26031d0b511e5d740dc3334a8c86923ac1b5 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Fri, 10 Apr 2026 11:31:40 +0100 Subject: [PATCH 15/23] Specify nla when printing ordering --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 561b5d1c7c..6603cccb3f 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -389,16 +389,16 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // no2hop improves the quality of ordering in general options[METIS_OPTION_NO2HOP] = 1; - logger_.printInfo("Running Metis\n"); + logger_.printInfo("Running Metis for %s\n", nla.c_str()); std::vector iperm(n); Int status = Highs_METIS_NodeND(&n, full_ptr.data(), full_rows.data(), NULL, options, permutations[i].data(), iperm.data()); - logger_.printInfo("Metis done\n"); + logger_.printInfo("Metis done for %s\n", nla.c_str()); if (status != METIS_OK) { - logger_.printInfo("Error with Metis\n"); + logger_.printInfo("Error with Metis for \n", nla.c_str()); failure[i] = true; } } else if (orderings_to_try[i] == kHipoAmdString) { @@ -406,23 +406,23 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, Highs_amd_defaults(control); double info[AMD_INFO]; - logger_.printInfo("Running AMD\n"); + logger_.printInfo("Running AMD for %s\n", nla.c_str()); Int status = Highs_amd_order(n, full_ptr.data(), full_rows.data(), permutations[i].data(), control, info); - logger_.printInfo("AMD done\n"); + logger_.printInfo("AMD done for %s\n", nla.c_str()); if (status != AMD_OK) { - logger_.printInfo("Error with AMD\n"); + logger_.printInfo("Error with AMD for %s\n", nla.c_str()); failure[i] = true; } } else if (orderings_to_try[i] == kHipoRcmString) { - logger_.printInfo("Running RCM\n"); + logger_.printInfo("Running RCM for %s\n", nla.c_str()); Int status = Highs_genrcm(n, full_ptr.back(), full_ptr.data(), full_rows.data(), permutations[i].data()); - logger_.printInfo("RCM done\n"); + logger_.printInfo("RCM done for %s\n", nla.c_str()); if (status != 0) { - logger_.printInfo("Error with RCM\n"); + logger_.printInfo("Error with RCM for %s\n", nla.c_str()); failure[i] = true; } } else { From 31f0badf3c6f22204f24e94ac23007ee10fd925e Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 13 Apr 2026 11:30:31 +0100 Subject: [PATCH 16/23] Fix warning --- highs/ipm/hipo/auxiliary/Logger.cpp | 2 -- highs/ipm/hipo/auxiliary/Logger.h | 3 +-- 2 files changed, 1 insertion(+), 4 deletions(-) diff --git a/highs/ipm/hipo/auxiliary/Logger.cpp b/highs/ipm/hipo/auxiliary/Logger.cpp index 3cc4688d44..2cbcade12f 100644 --- a/highs/ipm/hipo/auxiliary/Logger.cpp +++ b/highs/ipm/hipo/auxiliary/Logger.cpp @@ -2,8 +2,6 @@ namespace hipo { -Logger::Logger(Int level) : dev_level_{level} {} - void Logger::setOptions(const HighsLogOptions* log_options) { log_options_ = log_options; } diff --git a/highs/ipm/hipo/auxiliary/Logger.h b/highs/ipm/hipo/auxiliary/Logger.h index b670c90ee0..765d997699 100644 --- a/highs/ipm/hipo/auxiliary/Logger.h +++ b/highs/ipm/hipo/auxiliary/Logger.h @@ -12,11 +12,10 @@ namespace hipo { class Logger { - const Int dev_level_ = 0; const HighsLogOptions* log_options_; public: - Logger(Int level = 0); + Logger() = default; void setOptions(const HighsLogOptions* log_options); bool debug(Int level) const; From 9148f1ba9a82e9546d1e3b63a4f36765559abb09 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 13 Apr 2026 11:55:00 +0100 Subject: [PATCH 17/23] Fix overflow --- highs/ipm/hipo/ipm/Model.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/highs/ipm/hipo/ipm/Model.cpp b/highs/ipm/hipo/ipm/Model.cpp index 4cb73407a8..312bfd740c 100644 --- a/highs/ipm/hipo/ipm/Model.cpp +++ b/highs/ipm/hipo/ipm/Model.cpp @@ -48,8 +48,8 @@ void Model::nzBounds() { NE_nz_lb_ = A_.num_row_; NE_nz_ub_ = A_.num_row_; for (Int col = 0; col < A_.num_col_; ++col) { - Int used = 0; - Int unused = 0; + Int64 used = 0; + Int64 unused = 0; for (Int el = A_.start_[col]; el < A_.start_[col + 1]; ++el) { const Int row = A_.index_[el]; if (mark[row]) From 1f15d566c6193425f9bef2e631cb89b0378c9263 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 13 Apr 2026 11:55:18 +0100 Subject: [PATCH 18/23] Deal with failures and skips --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 32 +++++++++++++++++------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 6603cccb3f..9232449c9e 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -208,11 +208,17 @@ Int FactorHiGHSSolver::chooseNla() { bool overflow_NE = false; bool overflow_AS = false; - auto run_structure_NE = [&]() { - bool expect_AS_much_cheaper = - model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; + bool expect_AS_much_cheaper = + model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; + bool expect_NE_much_cheaper = + model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; + + bool can_skip_AS = !(model_.nonSeparableQp() || model_.m() == 0); + bool can_skip_NE = true; - if (expect_AS_much_cheaper || model_.nonSeparableQp() || model_.m() == 0) { + auto run_structure_NE = [&]() { + if ((expect_AS_much_cheaper && can_skip_NE) || model_.nonSeparableQp() || + model_.m() == 0) { failure_NE = true; logger_.printInfo("NE skipped\n"); } else { @@ -235,11 +241,6 @@ Int FactorHiGHSSolver::chooseNla() { }; auto run_analyse_AS = [&]() { - bool expect_NE_much_cheaper = - model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; - - bool can_skip_AS = !(model_.nonSeparableQp() || model_.m() == 0); - if (expect_NE_much_cheaper && can_skip_AS) { failure_AS = true; logger_.printInfo("AS skipped\n"); @@ -272,8 +273,21 @@ Int FactorHiGHSSolver::chooseNla() { tg.spawn([&]() { run_structure_NE(); }); tg.taskWait(); } + + // if NE was skipped but AS failed, use NE + if (expect_AS_much_cheaper && failure_AS) { + can_skip_NE = false; + run_structure_NE(); + } + run_analyse_NE(); + // if AS was skipped but NE failed, use AS + if (expect_NE_much_cheaper && failure_NE) { + can_skip_AS = false; + run_analyse_AS(); + } + Int status = kStatusOk; std::stringstream log_stream; From 1e95941fb3f809fbe04503de9caebee3ac750459 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 13 Apr 2026 15:16:28 +0100 Subject: [PATCH 19/23] Fix logic --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 27 ++++++++++++++---------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 9232449c9e..bfa8a81275 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -208,17 +208,22 @@ Int FactorHiGHSSolver::chooseNla() { bool overflow_NE = false; bool overflow_AS = false; - bool expect_AS_much_cheaper = + const bool AS_possible = true; + const bool NE_possible = !(model_.nonSeparableQp() || model_.m() == 0); + + const bool expect_AS_much_cheaper = model_.nzNElb() > model_.nzAS() * kNzBoundsRatio; - bool expect_NE_much_cheaper = + const bool expect_NE_much_cheaper = model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; - bool can_skip_AS = !(model_.nonSeparableQp() || model_.m() == 0); - bool can_skip_NE = true; + const bool skip_AS = NE_possible && expect_NE_much_cheaper; + const bool skip_NE = AS_possible && expect_AS_much_cheaper; + + bool allow_skip_AS = true; + bool allow_skip_NE = true; auto run_structure_NE = [&]() { - if ((expect_AS_much_cheaper && can_skip_NE) || model_.nonSeparableQp() || - model_.m() == 0) { + if (!NE_possible || (skip_NE && allow_skip_NE)) { failure_NE = true; logger_.printInfo("NE skipped\n"); } else { @@ -241,7 +246,7 @@ Int FactorHiGHSSolver::chooseNla() { }; auto run_analyse_AS = [&]() { - if (expect_NE_much_cheaper && can_skip_AS) { + if (!AS_possible || (skip_AS && allow_skip_AS)) { failure_AS = true; logger_.printInfo("AS skipped\n"); } else { @@ -275,16 +280,16 @@ Int FactorHiGHSSolver::chooseNla() { } // if NE was skipped but AS failed, use NE - if (expect_AS_much_cheaper && failure_AS) { - can_skip_NE = false; + if (skip_NE && failure_AS) { + allow_skip_NE = false; run_structure_NE(); } run_analyse_NE(); // if AS was skipped but NE failed, use AS - if (expect_NE_much_cheaper && failure_NE) { - can_skip_AS = false; + if (skip_AS && failure_NE) { + allow_skip_AS = false; run_analyse_AS(); } From bf686afd7222a0dcd4a6af6a760c84e7acfb6e3e Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Mon, 13 Apr 2026 15:18:44 +0100 Subject: [PATCH 20/23] Try multiple orderings --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 175 +++++++++++++---------- 1 file changed, 99 insertions(+), 76 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index bfa8a81275..5c9049bcf8 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -362,6 +362,10 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // - If ordering is "amd", "metis", "rcm" run only the ordering requested. // - If ordering is "choose", run "amd", "metis", and choose the best. + const std::string metis_2hop_string = "metis-2hop"; + const std::string metis_pfactor_string = "metis-pfactor"; + const std::string metis_2hop_pfactor_string = "metis-2hop-pfactor"; + // select which fill-reducing orderings should be tried std::vector orderings_to_try; if (options_.ordering != kHighsChooseString) @@ -369,6 +373,28 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, else { orderings_to_try.push_back(kHipoAmdString); orderings_to_try.push_back(kHipoMetisString); + + if (options_.parallel != kHighsOffString) { + Int available_threads = highs::parallel::num_threads(); + + // The first 3 threads are used to run NE structure, amd on AS, metis on + // AS. If there are more threads, run more orderings + available_threads -= 3; + + if (available_threads > 0) { + orderings_to_try.push_back(metis_2hop_string); + available_threads--; + } + if (available_threads > 0) { + orderings_to_try.push_back(metis_pfactor_string); + available_threads--; + } + if (available_threads > 0) { + orderings_to_try.push_back(metis_2hop_pfactor_string); + available_threads--; + } + } + // rcm is much worse in general, so no point in trying for now } @@ -394,61 +420,62 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, std::vector symbolics(orderings_to_try.size(), S); auto run_ordering_and_analyse = [&](Int i) { - // compute ordering - if (orderings_to_try[i] == kHipoMetisString) { + logger_.printInfo("Running %s for %s\n", orderings_to_try[i].c_str(), + nla.c_str()); + + if (orderings_to_try[i] == kHipoMetisString || + orderings_to_try[i] == metis_2hop_string || + orderings_to_try[i] == metis_pfactor_string || + orderings_to_try[i] == metis_2hop_pfactor_string) { idx_t options[METIS_NOPTIONS]; Highs_METIS_SetDefaultOptions(options); options[METIS_OPTION_SEED] = kMetisSeed; - // set logging of Metis depending on debug level options[METIS_OPTION_DBGLVL] = 0; - if (logger_.debug(2)) - options[METIS_OPTION_DBGLVL] = METIS_DBG_INFO | METIS_DBG_COARSEN; // no2hop improves the quality of ordering in general options[METIS_OPTION_NO2HOP] = 1; - logger_.printInfo("Running Metis for %s\n", nla.c_str()); - std::vector iperm(n); + if (orderings_to_try[i] == metis_2hop_string) + options[METIS_OPTION_NO2HOP] = 0; + else if (orderings_to_try[i] == metis_pfactor_string) + options[METIS_OPTION_PFACTOR] = 200; + else if (orderings_to_try[i] == metis_2hop_pfactor_string) { + options[METIS_OPTION_NO2HOP] = 0; + options[METIS_OPTION_PFACTOR] = 200; + } + std::vector iperm(n); Int status = Highs_METIS_NodeND(&n, full_ptr.data(), full_rows.data(), NULL, options, permutations[i].data(), iperm.data()); + if (status != METIS_OK) failure[i] = true; - logger_.printInfo("Metis done for %s\n", nla.c_str()); - if (status != METIS_OK) { - logger_.printInfo("Error with Metis for \n", nla.c_str()); - failure[i] = true; - } } else if (orderings_to_try[i] == kHipoAmdString) { double control[AMD_CONTROL]; Highs_amd_defaults(control); double info[AMD_INFO]; - - logger_.printInfo("Running AMD for %s\n", nla.c_str()); Int status = Highs_amd_order(n, full_ptr.data(), full_rows.data(), permutations[i].data(), control, info); - logger_.printInfo("AMD done for %s\n", nla.c_str()); + if (status != AMD_OK) failure[i] = true; - if (status != AMD_OK) { - logger_.printInfo("Error with AMD for %s\n", nla.c_str()); - failure[i] = true; - } } else if (orderings_to_try[i] == kHipoRcmString) { - logger_.printInfo("Running RCM for %s\n", nla.c_str()); Int status = Highs_genrcm(n, full_ptr.back(), full_ptr.data(), full_rows.data(), permutations[i].data()); - logger_.printInfo("RCM done for %s\n", nla.c_str()); + if (status != 0) failure[i] = true; - if (status != 0) { - logger_.printInfo("Error with RCM for %s\n", nla.c_str()); - failure[i] = true; - } } else { assert(1 == 0); } - if (failure[i]) return; + logger_.printInfo("Finished %s for %s\n", orderings_to_try[i].c_str(), + nla.c_str()); + + if (failure[i]) { + logger_.printInfo("Error with %s for %s\n", orderings_to_try[i].c_str(), + nla.c_str()); + return; + } failure[i] = FH_.analyse(symbolics[i], rows, ptr, signs, permutations[i]); @@ -471,63 +498,59 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, if (!b) ++num_success; } - if (orderings_to_try.size() < 2) { - S = std::move(symbolics[0]); - ordering = orderings_to_try[0]; - - } else if (orderings_to_try.size() == 2) { - // if there's only one success, obvious choice - if (failure[0] && !failure[1]) { - S = std::move(symbolics[1]); - ordering = orderings_to_try[1]; - } else if (!failure[0] && failure[1]) { - S = std::move(symbolics[0]); - ordering = orderings_to_try[0]; + if (num_success > 0) { + for (Int i = 0; i < static_cast(orderings_to_try.size()); ++i) { + if (!failure[i]) + logger_.printInfo( + "%20s for %s: %.2e %.2f\n", orderings_to_try[i].c_str(), + nla.c_str(), symbolics[i].flops(), + static_cast(symbolics[i].size()) / symbolics[i].sn()); + } + + // find the ordering with best flops + double best_flops = kHighsInf; + for (Int i = 0; i < static_cast(orderings_to_try.size()); ++i) { + if (!failure[i] && symbolics[i].flops() < best_flops) { + best_flops = symbolics[i].flops(); + } } - else if (num_success > 1) { - // need to choose the better ordering - - const double flops_0 = symbolics[0].flops(); - const double flops_1 = symbolics[1].flops(); - const double sn_avg_0 = symbolics[0].size() / symbolics[0].sn(); - const double sn_avg_1 = symbolics[1].size() / symbolics[1].sn(); - const double bytes_0 = symbolics[0].storage(); - const double bytes_1 = symbolics[1].storage(); - - Int chosen = -1; - - // selection rule: - // - if flops have a clear winner (+/- 20%), then choose it. - // - otherwise, choose the one with larger supernodes. - - if (flops_0 > kFlopsOrderingThresh * flops_1) - chosen = 1; - else if (flops_1 > kFlopsOrderingThresh * flops_0) - chosen = 0; - else if (sn_avg_0 > sn_avg_1) - chosen = 0; - else - chosen = 1; - - // fix selection if one or more require too much memory - const double bytes_thresh = kLargeStorageGB * 1024 * 1024 * 1024; - if (bytes_0 > bytes_thresh || bytes_1 > bytes_thresh) { - if (bytes_0 > bytes_1) - chosen = 1; - else - chosen = 0; + // find orderings with flops within kFlopsOrderingThresh of the best + std::vector consider; + for (Int i = 0; i < static_cast(orderings_to_try.size()); ++i) { + if (!failure[i] && + symbolics[i].flops() <= kFlopsOrderingThresh * best_flops) { + consider.push_back(i); } + } - assert(chosen == 0 || chosen == 1); + // among them, find the one with largest supernodes + double best_sn_size = 0.0; + Int chosen = -1; + for (Int i : consider) { + double sn_avg_size = + static_cast(symbolics[i].size()) / symbolics[i].sn(); + if (sn_avg_size > best_sn_size) { + best_sn_size = sn_avg_size; + chosen = i; + } + } - S = std::move(symbolics[chosen]); - ordering = orderings_to_try[chosen]; + // fix selection if one or more require too much memory + const double bytes_thresh = kLargeStorageGB * 1024 * 1024 * 1024; + double best_memory = kHighsInf; + Int ind_best_memory = -1; + for (Int i = 0; i < static_cast(orderings_to_try.size()); ++i) { + if (symbolics[i].storage() < best_memory) { + best_memory = symbolics[i].storage(); + ind_best_memory = i; + } } - } else { - // only two orderings tried for now - assert(0 == 1); + assert(chosen >= 0 && chosen < static_cast(orderings_to_try.size())); + + S = std::move(symbolics[chosen]); + ordering = orderings_to_try[chosen]; } return num_success > 0 ? kStatusOk : kStatusErrorAnalyse; From aef6d843a412ee2233224a0b1802ae146a528852 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Tue, 14 Apr 2026 12:22:54 +0100 Subject: [PATCH 21/23] Remove multiple orderings --- highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp | 39 +----------------------- 1 file changed, 1 insertion(+), 38 deletions(-) diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp index 5c9049bcf8..46da96bebb 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -362,10 +362,6 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // - If ordering is "amd", "metis", "rcm" run only the ordering requested. // - If ordering is "choose", run "amd", "metis", and choose the best. - const std::string metis_2hop_string = "metis-2hop"; - const std::string metis_pfactor_string = "metis-pfactor"; - const std::string metis_2hop_pfactor_string = "metis-2hop-pfactor"; - // select which fill-reducing orderings should be tried std::vector orderings_to_try; if (options_.ordering != kHighsChooseString) @@ -374,27 +370,6 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, orderings_to_try.push_back(kHipoAmdString); orderings_to_try.push_back(kHipoMetisString); - if (options_.parallel != kHighsOffString) { - Int available_threads = highs::parallel::num_threads(); - - // The first 3 threads are used to run NE structure, amd on AS, metis on - // AS. If there are more threads, run more orderings - available_threads -= 3; - - if (available_threads > 0) { - orderings_to_try.push_back(metis_2hop_string); - available_threads--; - } - if (available_threads > 0) { - orderings_to_try.push_back(metis_pfactor_string); - available_threads--; - } - if (available_threads > 0) { - orderings_to_try.push_back(metis_2hop_pfactor_string); - available_threads--; - } - } - // rcm is much worse in general, so no point in trying for now } @@ -423,10 +398,7 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, logger_.printInfo("Running %s for %s\n", orderings_to_try[i].c_str(), nla.c_str()); - if (orderings_to_try[i] == kHipoMetisString || - orderings_to_try[i] == metis_2hop_string || - orderings_to_try[i] == metis_pfactor_string || - orderings_to_try[i] == metis_2hop_pfactor_string) { + if (orderings_to_try[i] == kHipoMetisString) { idx_t options[METIS_NOPTIONS]; Highs_METIS_SetDefaultOptions(options); options[METIS_OPTION_SEED] = kMetisSeed; @@ -436,15 +408,6 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // no2hop improves the quality of ordering in general options[METIS_OPTION_NO2HOP] = 1; - if (orderings_to_try[i] == metis_2hop_string) - options[METIS_OPTION_NO2HOP] = 0; - else if (orderings_to_try[i] == metis_pfactor_string) - options[METIS_OPTION_PFACTOR] = 200; - else if (orderings_to_try[i] == metis_2hop_pfactor_string) { - options[METIS_OPTION_NO2HOP] = 0; - options[METIS_OPTION_PFACTOR] = 200; - } - std::vector iperm(n); Int status = Highs_METIS_NodeND(&n, full_ptr.data(), full_rows.data(), NULL, From a3f1367329715a90901081b7b92b2c7df96bb327 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Thu, 16 Apr 2026 15:19:00 +0100 Subject: [PATCH 22/23] Undo regularisation change to keep same behaviour --- highs/ipm/hipo/factorhighs/DenseFactKernel.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp index c30c44e2b1..eada1cdcf1 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp @@ -44,9 +44,9 @@ static void staticReg(double& pivot, Int sign, const Regul& regval, double old_pivot = pivot; if (sign > 0) - pivot += regval.dual; + pivot += regval.primal; else - pivot -= regval.primal; + pivot -= regval.dual; totalreg = pivot - old_pivot; } @@ -110,8 +110,8 @@ static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, assert(r >= 0); res = maxInCol(j, n, r, A, lda); double gamma_r = res.second; - double Arr = sign[r] > 0 ? A[r + lda * r] + regval.dual - : A[r + lda * r] - regval.primal; + double Arr = sign[r] > 0 ? A[r + lda * r] + regval.primal + : A[r + lda * r] - regval.dual; if ((std::abs(Ajj) >= kAlphaBK * gamma_j || std::abs(Ajj) * gamma_r >= kAlphaBK * gamma_j * gamma_j)) { @@ -198,7 +198,7 @@ Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, if (std::isnan(Ajj)) return kRetInvalidPivot; // add regularisation - //staticReg(Ajj, pivot_sign[j], regval, totalreg[j]); + staticReg(Ajj, pivot_sign[j], regval, totalreg[j]); #ifndef HIPO_PIVOTING // add static regularisation From 0688910dfa71fa6bbb99a5915124f5f70dcddab0 Mon Sep 17 00:00:00 2001 From: filikat <78378227+filikat@users.noreply.github.com> Date: Fri, 17 Apr 2026 11:42:47 +0100 Subject: [PATCH 23/23] Check pointer --- highs/ipm/hipo/ipm/Solver.cpp | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index 1f2400c62c..88bb5cb2f2 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -145,6 +145,10 @@ bool Solver::initialise() { // initialise linear solver LS_.reset(new FactorHiGHSSolver(*kkt_, options_, model_, regul_, info_, it_->data, logger_)); + if (!LS_) { + info_.status = kStatusError; + return true; + } if (Int status = LS_->setup()) { info_.status = (Status)status; return true;