diff --git a/cmake/sources.cmake b/cmake/sources.cmake index 89897224c9..d789272aef 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -182,7 +182,7 @@ set(hipo_sources ipm/hipo/ipm/FactorHiGHSSolver.cpp ipm/hipo/ipm/Control.cpp ipm/hipo/ipm/Iterate.cpp - ipm/hipo/ipm/LogHighs.cpp + ipm/hipo/ipm/KktMatrix.cpp ipm/hipo/ipm/Model.cpp ipm/hipo/ipm/PreProcess.cpp ipm/hipo/ipm/Refine.cpp @@ -195,8 +195,8 @@ 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 ipm/hipo/ipm/Options.h ipm/hipo/ipm/PreProcess.h @@ -246,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/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/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..2cbcade12f --- /dev/null +++ b/highs/ipm/hipo/auxiliary/Logger.cpp @@ -0,0 +1,42 @@ +#include "Logger.h" + +namespace hipo { + +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..765d997699 --- /dev/null +++ b/highs/ipm/hipo/auxiliary/Logger.h @@ -0,0 +1,79 @@ +#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 HighsLogOptions* log_options_; + + public: + Logger() = default; + 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 587c71a845..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. @@ -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() { @@ -941,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"); } } } @@ -1341,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 9a8936421a..eada1cdcf1 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 af06414b82..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; @@ -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) { @@ -197,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 @@ -248,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; } @@ -326,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 b0ba1d7682..46da96bebb 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp @@ -5,29 +5,25 @@ #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" 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), + 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}, - 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() { @@ -35,225 +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; - - ptrAS_.resize(nA_ + mA_ + 1); - rowsAS_.resize(nzBlock11 + nzA_ + mA_); - valAS_.resize(nzBlock11 + nzA_ + mA_); - - Int next = 0; - - for (Int i = 0; i < nA_; ++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] = nA_ + 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 < mA_; ++i) { - rowsAS_[next] = nA_ + i; - ++next; - ptrAS_[nA_ + i + 1] = 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(!ptrAS_.empty() && !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); - } - - 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]++; - } - } - } - - ptrNE_.clear(); - rowsNE_.clear(); - - // ptr is allocated its exact size - 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)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; - } - } - - 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(!ptrNE_.empty() && !rowsNE_.empty()); - - valNE_.resize(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 = ptrNE_[row]; el < ptrNE_[row + 1]; ++el) { - Int index = rowsNE_[el]; - valNE_[el] = work[index]; - work[index] = 0.0; - } - } - - return kStatusOk; -} - // ========================================================================= // Analyse phase // ========================================================================= @@ -262,17 +39,20 @@ 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; + + 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"); + logger_.printInfo("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,16 +63,16 @@ 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); + 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(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; @@ -306,18 +86,14 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { // only execute factorisation if it has not been done yet assert(!this->valid_); - Clock clock; - - // build matrix - buildASvalues(scaling); - info_.matrix_time += clock.stop(); + kkt_.buildASvalues(scaling); // set static regularisation, since it may have changed FH_.setRegularisation(regul_.primal, regul_.dual); - // factorise matrix - clock.start(); - if (FH_.factorise(S_, rowsAS_, ptrAS_, valAS_)) return kStatusErrorFactorise; + Clock clock; + if (FH_.factorise(S_, kkt_.rowsAS, kkt_.ptrAS, kkt_.valAS)) + return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -329,19 +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 - buildNEvalues(scaling); - info_.matrix_time += clock.stop(); + kkt_.buildNEvalues(scaling); // set static regularisation, since it may have changed FH_.setRegularisation(regul_.primal, regul_.dual); - // factorise - clock.start(); - if (FH_.factorise(S_, rowsNE_, ptrNE_, valNE_)) return kStatusErrorFactorise; - + Clock clock; + if (FH_.factorise(S_, kkt_.rowsNE, kkt_.ptrNE, kkt_.valNE)) + return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -406,17 +177,23 @@ Int FactorHiGHSSolver::solveNE(const std::vector& rhs, // ========================================================================= Int FactorHiGHSSolver::setup() { + Clock clock; + if (Int status = setNla()) return status; setParallel(); - S_.print(log_, log_.debug(1)); + std::stringstream log_stream; + log_stream << textline("Analyse time:") << fix(clock.stop(), 0, 2) << '\n'; + logger_.print(log_stream.str().c_str()); + + 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; } @@ -431,17 +208,30 @@ Int FactorHiGHSSolver::chooseNla() { bool overflow_NE = false; bool overflow_AS = false; + 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; + const bool expect_NE_much_cheaper = + model_.nzAS() > model_.nzNEub() * kNzBoundsRatio; + + 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 ((model_.m() > kMinRowsForDensity && - model_.maxColDensity() > kDenseColThresh) || - model_.nonSeparableQp() || model_.m() == 0) { + if (!NE_possible || (skip_NE && allow_skip_NE)) { failure_NE = true; + logger_.printInfo("NE skipped\n"); } else { - Int status = buildNEstructure(); + 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; @@ -456,19 +246,24 @@ Int FactorHiGHSSolver::chooseNla() { }; auto run_analyse_AS = [&]() { - Int AS_status = 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 (!AS_possible || (skip_AS && allow_skip_AS)) { + failure_AS = true; + 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) { + logger_.printInfo("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; - 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 @@ -483,8 +278,21 @@ Int FactorHiGHSSolver::chooseNla() { tg.spawn([&]() { run_structure_NE(); }); tg.taskWait(); } + + // if NE was skipped but AS failed, use NE + 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 (skip_AS && failure_NE) { + allow_skip_AS = false; + run_analyse_AS(); + } + Int status = kStatusOk; std::stringstream log_stream; @@ -502,7 +310,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 @@ -530,15 +338,15 @@ Int FactorHiGHSSolver::chooseNla() { } } - log_.print(log_stream); + logger_.print(log_stream.str().c_str()); if (status == kStatusOk) { if (options_.nla == kHipoAugmentedString) { S_ = std::move(symb_AS); - freeNEmemory(); + kkt_.freeNEmemory(); } else { S_ = std::move(symb_NE); - freeASmemory(); + kkt_.freeASmemory(); } } @@ -554,15 +362,14 @@ 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) 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 } @@ -570,8 +377,8 @@ 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)) { - log_.printDevInfo("NE interrupted\n"); + if (ptr.back() >= kkt_.NE_nz_limit.load(std::memory_order_relaxed)) { + logger_.printInfo("NE interrupted\n"); return kStatusErrorAnalyse; } } @@ -588,67 +395,56 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, std::vector symbolics(orderings_to_try.size(), S); auto run_ordering_and_analyse = [&](Int i) { - // compute ordering + logger_.printInfo("Running %s for %s\n", orderings_to_try[i].c_str(), + nla.c_str()); + if (orderings_to_try[i] == kHipoMetisString) { 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 (log_.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"); 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; - log_.printDevInfo("Metis done\n"); - if (status != METIS_OK) { - log_.printDevInfo("Error with Metis\n"); - failure[i] = true; - } } else if (orderings_to_try[i] == kHipoAmdString) { double control[AMD_CONTROL]; Highs_amd_defaults(control); double info[AMD_INFO]; - - log_.printDevInfo("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"); + if (status != AMD_OK) failure[i] = true; - if (status != AMD_OK) { - log_.printDevInfo("Error with AMD\n"); - failure[i] = true; - } } else if (orderings_to_try[i] == kHipoRcmString) { - log_.printDevInfo("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"); + if (status != 0) failure[i] = true; - if (status != 0) { - log_.printDevInfo("Error with RCM\n"); - 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]); - 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); } }; @@ -665,63 +461,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()); } - 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 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(); } + } - assert(chosen == 0 || chosen == 1); + // 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); + } + } - S = std::move(symbolics[chosen]); - ordering = orderings_to_try[chosen]; + // 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; + } } - } else { - // only two orderings tried for now - assert(0 == 1); + // 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; + } + } + + 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; @@ -731,30 +523,30 @@ 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; } 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"); + 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"; } 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"); + 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"; @@ -765,12 +557,14 @@ 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(); return kStatusOk; } @@ -852,7 +646,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); } @@ -867,21 +661,4 @@ void FactorHiGHSSolver::getReg(std::vector& reg) { FH_.getRegularisation(reg); } -void FactorHiGHSSolver::freeASmemory() { - // Give up memory used for AS. - freeVector(ptrAS_); - freeVector(rowsAS_); - freeVector(valAS_); -} - -void FactorHiGHSSolver::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/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h index 4e2d88328b..651147aadb 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" @@ -20,31 +21,15 @@ class FactorHiGHSSolver : public LinearSolver { // symbolic factorisation Symbolic S_; - // 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_; + KktMatrix& kkt_; const Regularisation& regul_; - Info& info_; IpmData& data_; - const LogHighs& log_; - + const Logger& logger_; 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"; @@ -54,22 +39,13 @@ 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); 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); + const Logger& logger); // Override functions Int factorAS(const std::vector& scaling) override; 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/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_) { diff --git a/highs/ipm/hipo/ipm/KktMatrix.cpp b/highs/ipm/hipo/ipm/KktMatrix.cpp new file mode 100644 index 0000000000..304ea874ff --- /dev/null +++ b/highs/ipm/hipo/ipm/KktMatrix.cpp @@ -0,0 +1,265 @@ +#include "KktMatrix.h" + +#include "ipm/hipo/auxiliary/Auxiliary.h" + +namespace hipo { + +KktMatrix::KktMatrix(const Model& m, const Regularisation& r, Info& i, + const Logger& l) + : model{m}, regul{r}, info{i}, logger{l} {} + +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(); + const Int n = A.num_col_; + const Int m = A.num_row_; + const Int nzA = A.numNz(); + const Int nzQ = Q.numNz(); + + logger.printInfo("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; + } + + info.AS_structure_time = clock.stop(); + + return kStatusOk; +} + +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_; + + 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); + } + + info.matrix_time += clock.stop(); + + 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 + Clock clock; + + const HighsSparseMatrix& A = model.A(); + const Int n = A.num_col_; + const Int m = A.num_row_; + const Int nzA = A.numNz(); + + logger.printInfo("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; + } + } + + info.NE_structure_time = clock.stop(); + 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()); + Clock clock; + + 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; + } + } + + info.matrix_time += clock.stop(); + + 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..c83ea75402 --- /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 Logger& logger; + + KktMatrix(const Model& model, const Regularisation& regul, Info& info, + const Logger& logger); + + 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/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 c1a6fbf6d1..312bfd740c 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 { @@ -27,7 +27,7 @@ Int Model::init(const HighsLp& lp, const HighsHessian& Q) { preprocess(); - denseColumns(); + nzBounds(); computeNorms(); // double transpose to sort indices of each column @@ -41,6 +41,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) { + 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]) + ++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. @@ -150,15 +177,13 @@ 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'; 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()) @@ -252,12 +277,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; @@ -280,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"; @@ -315,24 +343,17 @@ 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); - - log.print(log_stream); -} - -void Model::denseColumns() { - // Compute the maximum density of any column of A and count the number of - // dense columns. + if (logger.debug(1)) { + preprocessor_.print(log_stream); + log_stream << textline("Free variables:") << integer(free_vars) << '\n'; - 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_; + 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"; } + + 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 11dec6e4bd..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" @@ -49,8 +49,8 @@ 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_; @@ -66,7 +66,6 @@ class Model { inf_norm_rows_; void preprocess(); - void denseColumns(); Int checkData() const; void computeNorms(); @@ -75,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; @@ -89,6 +88,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]); } @@ -115,12 +116,13 @@ 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]; } 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; diff --git a/highs/ipm/hipo/ipm/Parameters.h b/highs/ipm/hipo/ipm/Parameters.h index b82e10494d..9ad471b313 100644 --- a/highs/ipm/hipo/ipm/Parameters.h +++ b/highs/ipm/hipo/ipm/Parameters.h @@ -39,9 +39,8 @@ 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 const Int kMaxIterRefine = 3; diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index f159e042e1..88bb5cb2f2 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) { @@ -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,9 +136,19 @@ bool Solver::initialise() { start_time_ = control_.elapsed(); + kkt_.reset(new KktMatrix(model_, regul_, info_, logger_)); + if (!kkt_) { + info_.status = kStatusError; + return true; + } + // 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, logger_)); + if (!LS_) { + info_.status = kStatusError; + return true; + } if (Int status = LS_->setup()) { info_.status = (Status)status; return true; @@ -164,9 +178,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() { @@ -252,7 +263,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; } @@ -264,7 +275,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; } @@ -275,9 +286,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; } @@ -299,7 +310,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; @@ -312,11 +323,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; } @@ -336,7 +347,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; } @@ -345,7 +356,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; } @@ -370,7 +381,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; } @@ -379,7 +390,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; } @@ -606,13 +617,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; } @@ -623,7 +634,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; } @@ -632,7 +643,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; } @@ -691,7 +702,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; } @@ -705,7 +716,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; } @@ -985,11 +996,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; } @@ -1000,7 +1011,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; } } @@ -1011,7 +1022,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; } @@ -1039,20 +1050,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; @@ -1073,12 +1084,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; } @@ -1086,7 +1097,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; } } @@ -1098,14 +1109,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; @@ -1118,10 +1129,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; } @@ -1138,58 +1149,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 { @@ -1203,10 +1193,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 { @@ -1223,7 +1213,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'; @@ -1255,7 +1245,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 7254a957e2..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" @@ -32,6 +32,8 @@ class Solver { // Linear solver interface std::unique_ptr LS_; + std::unique_ptr kkt_; + // Iterate object interface std::unique_ptr it_; @@ -64,7 +66,7 @@ class Solver { ipx::LpSolver ipx_lps_; // Interface to Highs logging - LogHighs logH_; + Logger logger_; double start_time_;