diff --git a/.github/workflows/hipo-linear-solver.yml b/.github/workflows/hipo-linear-solver.yml new file mode 100644 index 00000000000..0e35c03c2ca --- /dev/null +++ b/.github/workflows/hipo-linear-solver.yml @@ -0,0 +1,37 @@ +name: hipo-linear-solver + +on: [push, pull_request] + +jobs: + ubuntu: + runs-on: ubuntu-latest + strategy: + fail-fast: false + + steps: + - uses: actions/checkout@v6 + + - name: Create Build Environment + run: cmake -E make_directory ${{github.workspace}}/build + + - name: Configure CMake + working-directory: ${{github.workspace}}/build + run: | + cmake $GITHUB_WORKSPACE -DHIPO=ON -DBUILD_OPENBLAS=ON + + - name: Build + working-directory: ${{github.workspace}}/build + run: | + cmake --build . -j2 + + - name: Compile example + working-directory: ${{github.workspace}} + run: | + clang highs/ipm/hipo/factorhighs/FactorHighs_c_api_example.c \ + -o example -Ihighs -Ibuild -Lbuild/lib -lhighs + + - name: Run example + working-directory: ${{github.workspace}} + run: | + export LD_LIBRARY_PATH=build/lib + ./example \ No newline at end of file diff --git a/check/TestHipo.cpp b/check/TestHipo.cpp index 340a6acc07e..389804b6615 100644 --- a/check/TestHipo.cpp +++ b/check/TestHipo.cpp @@ -6,6 +6,7 @@ #include "Highs.h" #include "catch.hpp" #include "io/Filereader.h" +#include "ipm/hipo/factorhighs/FactorHighs_c_api.h" #include "ipm/hipo/ipm/Solver.h" #include "ipm/hipo/ipm/Status.h" #include "lp_data/HighsCallback.h" @@ -127,4 +128,81 @@ TEST_CASE("test-hipo-freevar", "[highs_hipo]") { runHipoTest(highs, "perold.mps", -9.381e3, expected_status, "on"); runHipoTest(highs, "perold.mps", -9.381e3, expected_status, "off"); highs.resetGlobalScheduler(true); +} + +TEST_CASE("test-hipo-linear-solver", "[highs_hipo]") { + // problem size + const int n = 5; + const int nz = 10; + + // define the lower triangle in CSC format, with 0-based indexing + int ptr[n + 1] = {0, 3, 5, 8, 9, 10}; + int rows[nz] = {0, 2, 3, 1, 2, 2, 3, 4, 3, 4}; + double vals[nz] = {5, 3, 4, 3, 2, 9, -1, 1, 8, 1}; + + // matrix is spd, so expect all pivots to be positive + int signs[n] = {1, 1, 1, 1, 1}; + + // rhs and expected solution + double rhs[n] = {1, 2, 3, 4, 5}; + const double lhs[n] = {0.457627118644068, 1.118644067796610, + -0.677966101694915, 0.186440677966102, + 5.677966101694915}; + + // initialise with default number of threads + const int num_threads = 0; + int initialise_status = FactorHighs_initialise(num_threads); + if (initialise_status != 0) return; + + void* S = FactorHighs_symbolic_create(); + void* FH = FactorHighs_create(); + + // set options + const int logging_on = dev_run; + FactorHighs_setLogging(FH, logging_on); + + const int block_size = 64; + FactorHighs_setBlockSize(FH, block_size); + + const int pivoting_off = 0; + FactorHighs_setPivoting(FH, pivoting_off); + + FactorHighs_setRegularisation(FH, 0.0, 0.0); + + // compute ordering with metis + int perm[n]; + int metis_status = FactorHighs_reorderMetis(FH, n, nz, rows, ptr, perm); + REQUIRE(metis_status == 0); + + // perform analyse phase + int analyse_status = + FactorHighs_analyse(FH, S, n, nz, rows, ptr, signs, perm); + REQUIRE(analyse_status == 0); + + // print extended statistics of symbolic factorisation + int verbose = 1; + FactorHighs_symbolic_print(FH, S, verbose); + + // factorise the matrix + int factorise_status = FactorHighs_factorise(FH, S, n, nz, rows, ptr, vals); + REQUIRE(factorise_status == 0); + + // compute the inertia of the factorisation + int pos, neg, zero; + double zero_tolerance = 1e-16; + FactorHighs_inertia(FH, &pos, &neg, &zero, zero_tolerance); + + // triangular solve + int solve_status = FactorHighs_solve(FH, rhs, 1); + REQUIRE(solve_status == 0); + + // compute error + double error = 0.0; + for (int i = 0; i < n; ++i) error += fabs(rhs[i] - lhs[i]); + REQUIRE(error < 1e-12); + + // terminate + FactorHighs_symbolic_destroy(S); + FactorHighs_destroy(FH); + FactorHighs_terminate(); } \ No newline at end of file diff --git a/cmake/sources.cmake b/cmake/sources.cmake index a1683a8315a..1f558fefe57 100644 --- a/cmake/sources.cmake +++ b/cmake/sources.cmake @@ -260,7 +260,7 @@ set(ipx_headers set(hipo_sources ipm/hipo/ipm/IpmData.cpp - ipm/hipo/ipm/FactorHiGHSSolver.cpp + ipm/hipo/ipm/FactorHighsSolver.cpp ipm/hipo/ipm/Control.cpp ipm/hipo/ipm/Iterate.cpp ipm/hipo/ipm/KktMatrix.cpp @@ -271,7 +271,7 @@ set(hipo_sources set(hipo_headers ipm/hipo/ipm/IpmData.h - ipm/hipo/ipm/FactorHiGHSSolver.h + ipm/hipo/ipm/FactorHighsSolver.h ipm/hipo/ipm/Parameters.h ipm/hipo/ipm/Control.h ipm/hipo/ipm/Info.h @@ -292,7 +292,8 @@ set(factor_highs_sources ipm/hipo/factorhighs/DenseFactHybrid.cpp ipm/hipo/factorhighs/DenseFactKernel.cpp ipm/hipo/factorhighs/DgemmParallel.cpp - ipm/hipo/factorhighs/FactorHiGHS.cpp + ipm/hipo/factorhighs/FactorHighs_c_api.cpp + ipm/hipo/factorhighs/FactorHighs.cpp ipm/hipo/factorhighs/Factorise.cpp ipm/hipo/factorhighs/FormatHandler.cpp ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp @@ -310,8 +311,10 @@ set(factor_highs_headers ipm/hipo/factorhighs/DataCollector.h ipm/hipo/factorhighs/DenseFact.h ipm/hipo/factorhighs/DgemmParallel.h - ipm/hipo/factorhighs/FactorHiGHS.h - ipm/hipo/factorhighs/FactorHiGHSSettings.h + ipm/hipo/factorhighs/FactorHighs_c_api.h + ipm/hipo/factorhighs/FactorHighs.h + ipm/hipo/factorhighs/FactorHighsOptions.h + ipm/hipo/factorhighs/FactorHighsSettings.h ipm/hipo/factorhighs/Factorise.h ipm/hipo/factorhighs/FormatHandler.h ipm/hipo/factorhighs/HybridHybridFormatHandler.h diff --git a/extern/blas/mycblas.h b/extern/blas/mycblas.h index 85fd4a3506a..adc9712cea0 100644 --- a/extern/blas/mycblas.h +++ b/extern/blas/mycblas.h @@ -82,7 +82,6 @@ void openblas_set_num_threads(int num_threads); #endif void highs_openblas_set_num_threads(int num_threads); - #ifdef __cplusplus } #endif diff --git a/highs/ipm/hipo/auxiliary/Auxiliary.cpp b/highs/ipm/hipo/auxiliary/Auxiliary.cpp index 3d2a004e989..3243ae04594 100644 --- a/highs/ipm/hipo/auxiliary/Auxiliary.cpp +++ b/highs/ipm/hipo/auxiliary/Auxiliary.cpp @@ -26,21 +26,19 @@ void subtreeSize(const std::vector& parent, std::vector& sizes) { } } -void transpose(const std::vector& ptr, const std::vector& rows, - std::vector& ptrT, std::vector& rowsT) { +void transpose(Int n, Int nz, const Int* ptr, const Int* rows, Int* ptrT, + Int* rowsT) { // Compute the transpose of the matrix and return it in rowsT and ptrT - Int n = ptr.size() - 1; - std::vector work(n); // count the entries in each row into work - for (Int i = 0; i < ptr.back(); ++i) { + for (Int i = 0; i < nz; ++i) { ++work[rows[i]]; } // sum row sums to obtain pointers - counts2Ptr(ptrT, work); + counts2Ptr(n, ptrT, work.data()); for (Int j = 0; j < n; ++j) { for (Int el = ptr[j]; el < ptr[j + 1]; ++el) { @@ -54,21 +52,24 @@ 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) { - // Compute the transpose of the matrix and return it in rowsT, ptrT and valT + std::vector& ptrT, std::vector& rowsT) { + transpose(ptr.size() - 1, rows.size(), ptr.data(), rows.data(), ptrT.data(), + rowsT.data()); +} - Int n = ptr.size() - 1; +void transpose(Int n, Int nz, const Int* ptr, const Int* rows, + const double* val, Int* ptrT, Int* rowsT, double* valT) { + // Compute the transpose of the matrix and return it in rowsT, ptrT and valT std::vector work(n); // count the entries in each row into work - for (Int i = 0; i < ptr.back(); ++i) { + for (Int i = 0; i < nz; ++i) { ++work[rows[i]]; } // sum row sums to obtain pointers - counts2Ptr(ptrT, work); + counts2Ptr(n, ptrT, work.data()); for (Int j = 0; j < n; ++j) { for (Int el = ptr[j]; el < ptr[j + 1]; ++el) { @@ -82,11 +83,18 @@ 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) { + transpose(ptr.size() - 1, rows.size(), ptr.data(), rows.data(), val.data(), + ptrT.data(), rowsT.data(), valT.data()); +} + 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. + // 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); @@ -117,7 +125,7 @@ void permuteSym(const std::vector& iperm, std::vector& ptr, // get column pointers by summing the count of nonzeros in each column. // copy column pointers into work - counts2Ptr(new_ptr, work); + counts2Ptr(n, new_ptr.data(), work.data()); std::vector new_rows(new_ptr.back()); std::vector new_val; @@ -300,16 +308,15 @@ Int maxDepthTree(const std::vector& parent) { return max_depth; } -void fullFromLower(const std::vector& ptrL, const std::vector& rowsL, +void fullFromLower(Int n, Int nz, const Int* ptrL, const Int* rowsL, std::vector& ptrF, std::vector& rowsF) { // Given a sparse matrix in lower triangular format, build the same matrix in // full format, without diagonal entries. - std::vector rowsU(rowsL.size()); - std::vector ptrU(ptrL.size()); - transpose(ptrL, rowsL, ptrU, rowsU); + std::vector rowsU(nz); + std::vector ptrU(n + 1); + transpose(n, nz, ptrL, rowsL, ptrU.data(), rowsU.data()); - const Int n = ptrL.size() - 1; std::vector work(n); for (Int j = 0; j < n; ++j) { for (Int el = ptrU[j]; el < ptrU[j + 1]; ++el) { @@ -321,7 +328,7 @@ void fullFromLower(const std::vector& ptrL, const std::vector& rowsL, } ptrF.assign(n + 1, 0); - counts2Ptr(ptrF, work); + counts2Ptr(n, ptrF.data(), work.data()); rowsF.assign(ptrF.back(), 0); for (Int j = 0; j < n; ++j) { for (Int el = ptrU[j]; el < ptrU[j + 1]; ++el) { diff --git a/highs/ipm/hipo/auxiliary/Auxiliary.h b/highs/ipm/hipo/auxiliary/Auxiliary.h index 2a7b94f5f7f..b41d0b02f04 100644 --- a/highs/ipm/hipo/auxiliary/Auxiliary.h +++ b/highs/ipm/hipo/auxiliary/Auxiliary.h @@ -14,8 +14,12 @@ namespace hipo { void inversePerm(const std::vector& perm, std::vector& iperm); void subtreeSize(const std::vector& parent, std::vector& sizes); +void transpose(Int n, Int nz, const Int* ptr, const Int* rows, Int* ptrT, + Int* rowsT); void transpose(const std::vector& ptr, const std::vector& rows, std::vector& ptrT, std::vector& rowsT); +void transpose(Int n, Int nz, const Int* ptr, const Int* rows, + const double* val, Int* ptrT, Int* rowsT, double* valT); void transpose(const std::vector& ptr, const std::vector& rows, const std::vector& val, std::vector& ptrT, std::vector& rowsT, std::vector& valT); @@ -32,19 +36,18 @@ void processEdge(Int j, Int i, const std::vector& first, Int64 getDiagStart(Int n, Int k, Int nb, Int n_blocks, std::vector& start, bool triang = false); Int maxDepthTree(const std::vector& parent); -void fullFromLower(const std::vector& ptrL, const std::vector& rowsL, +void fullFromLower(Int n, Int nz, const Int* ptrL, const Int* rowsL, std::vector& ptrF, std::vector& rowsF); double snFlops(double size, double clique_size); double snSpops(double clique_size); template -void counts2Ptr(std::vector& ptr, std::vector& w) { +void counts2Ptr(Int n, T* ptr, T* w) { // Given the column counts in the vector w (of size n), // compute the column pointers in the vector ptr (of size n+1), // and copy the first n pointers back into w. T temp_nz{}; - T n = w.size(); for (T j = 0; j < n; ++j) { ptr[j] = temp_nz; temp_nz += w[j]; @@ -53,6 +56,14 @@ void counts2Ptr(std::vector& ptr, std::vector& w) { ptr[n] = temp_nz; } +template +void permuteVector(T* v, const std::vector& perm) { + // Permute vector v according to permutation perm. + const Int n = perm.size(); + std::vector temp_v(v, v + n); + for (Int i = 0; i < n; ++i) v[i] = temp_v[perm[i]]; +} + template void permuteVector(std::vector& v, const std::vector& perm) { // Permute vector v according to permutation perm. @@ -60,6 +71,14 @@ void permuteVector(std::vector& v, const std::vector& perm) { for (Int i = 0; i < static_cast(v.size()); ++i) v[i] = temp_v[perm[i]]; } +template +void permuteVectorInverse(T* v, const std::vector& iperm) { + // Permute vector v according to inverse permutation iperm. + const Int n = iperm.size(); + std::vector temp_v(v, v + n); + for (Int i = 0; i < n; ++i) v[iperm[i]] = temp_v[i]; +} + template void permuteVectorInverse(std::vector& v, const std::vector& iperm) { // Permute vector v according to inverse permutation iperm. diff --git a/highs/ipm/hipo/auxiliary/Logger.cpp b/highs/ipm/hipo/auxiliary/Logger.cpp index 2cbcade12fd..6eda9463fb6 100644 --- a/highs/ipm/hipo/auxiliary/Logger.cpp +++ b/highs/ipm/hipo/auxiliary/Logger.cpp @@ -1,7 +1,11 @@ #include "Logger.h" +#include + namespace hipo { +Logger::Logger(bool use_printf) : use_printf_{use_printf} {} + void Logger::setOptions(const HighsLogOptions* log_options) { log_options_ = log_options; } @@ -39,4 +43,11 @@ std::string integer(Int i, Int width) { return s.str(); } +void logger_printf(const char* format, ...) { + va_list argptr; + va_start(argptr, format); + vprintf(format, argptr); + va_end(argptr); +} + } // namespace hipo diff --git a/highs/ipm/hipo/auxiliary/Logger.h b/highs/ipm/hipo/auxiliary/Logger.h index 765d997699f..86516485288 100644 --- a/highs/ipm/hipo/auxiliary/Logger.h +++ b/highs/ipm/hipo/auxiliary/Logger.h @@ -11,11 +11,14 @@ namespace hipo { +void logger_printf(const char* format, ...); + class Logger { - const HighsLogOptions* log_options_; + const HighsLogOptions* log_options_ = nullptr; + const bool use_printf_ = false; public: - Logger() = default; + Logger(bool use_printf = false); void setOptions(const HighsLogOptions* log_options); bool debug(Int level) const; @@ -23,36 +26,48 @@ class Logger { void print(const char* format, Args... args) const { if (log_options_) highsLogUser(*log_options_, HighsLogType::kInfo, format, args...); + else if (use_printf_) + logger_printf(format, args...); } template void printw(const char* format, Args... args) const { if (log_options_) highsLogUser(*log_options_, HighsLogType::kWarning, format, args...); + else if (use_printf_) + logger_printf(format, args...); } template void printe(const char* format, Args... args) const { if (log_options_) highsLogUser(*log_options_, HighsLogType::kError, format, args...); + else if (use_printf_) + logger_printf(format, args...); } template void printInfo(const char* format, Args... args) const { if (log_options_) highsLogDev(*log_options_, HighsLogType::kInfo, format, args...); + else if (use_printf_) + logger_printf(format, args...); } template void printDetailed(const char* format, Args... args) const { if (log_options_) highsLogDev(*log_options_, HighsLogType::kDetailed, format, args...); + else if (use_printf_) + logger_printf(format, args...); } template void printVerbose(const char* format, Args... args) const { if (log_options_) highsLogDev(*log_options_, HighsLogType::kVerbose, format, args...); + else if (use_printf_) + logger_printf(format, args...); } }; diff --git a/highs/ipm/hipo/factorhighs/Analyse.cpp b/highs/ipm/hipo/factorhighs/Analyse.cpp index be474c81b3f..44531603a2e 100644 --- a/highs/ipm/hipo/factorhighs/Analyse.cpp +++ b/highs/ipm/hipo/factorhighs/Analyse.cpp @@ -8,7 +8,7 @@ #include #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "ReturnValues.h" #include "ipm/hipo/auxiliary/Auxiliary.h" #include "ipm/hipo/auxiliary/Logger.h" @@ -17,25 +17,36 @@ namespace hipo { 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 Logger* logger, - DataCollector& data, const std::vector& perm) - : logger_{logger}, data_{data} { +Analyse::Analyse(Int n, Int nz, const Int* rows, const Int* ptr, + const Int* signs, const FHoptions& FH_opt, + const Logger* logger, DataCollector& data, const Int* perm) + : FH_opt_{FH_opt}, 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. // Only the lower triangular part is used. + // Diagonal entries must be stored for each column, even if zero. // signs contains the sign that each pivot should have. - n_ = ptr.size() - 1; - nz_ = rows.size(); - signs_ = signs; - nb_ = nb; + n_ = n; + nz_ = nz; + + rows_lower_ = std::vector(rows, rows + nz_); + ptr_lower_ = std::vector(ptr, ptr + n_ + 1); + signs_ = std::vector(signs, signs + n_); + perm_ = std::vector(perm, perm + n_); + + // adjust data if one-based indexing is used + if (FH_opt_.one_indexing) { + for (Int& i : rows_lower_) --i; + for (Int& i : ptr_lower_) --i; + for (Int& i : perm_) --i; + } // Create upper triangular part rows_upper_.resize(nz_); ptr_upper_.resize(n_ + 1); - transpose(ptr, rows, ptr_upper_, rows_upper_); + transpose(ptr_lower_, rows_lower_, ptr_upper_, rows_upper_); // Permute the matrix with identical permutation, to extract upper triangular // part, if the input is not lower triangular. @@ -55,7 +66,6 @@ Analyse::Analyse(const std::vector& rows, const std::vector& ptr, transpose(ptr_upper_, rows_upper_, ptr_lower_, rows_lower_); transpose(ptr_lower_, rows_lower_, ptr_upper_, rows_upper_); - perm_ = perm; iperm_.resize(n_); inversePerm(perm_, iperm_); @@ -724,7 +734,7 @@ void Analyse::snPattern() { std::vector work(sn_indices_.size()); for (Int i = 0; i < static_cast(sn_indices_.size()); ++i) work[i] = sn_indices_[i]; - counts2Ptr(ptr_sn_, work); + counts2Ptr(work.size(), ptr_sn_.data(), work.data()); // consider each row for (Int i = 0; i < n_; ++i) { @@ -889,17 +899,18 @@ void Analyse::computeStorage(Int fr, Int sz, Int64& fr_entries, // compute storage required by frontal and clique, based on the format used const Int cl = fr - sz; + const Int nb = FH_opt_.nb; - Int n_blocks = (sz - 1) / nb_ + 1; + Int n_blocks = (sz - 1) / nb + 1; std::vector temp; - fr_entries = getDiagStart(fr, sz, nb_, n_blocks, temp); + fr_entries = getDiagStart(fr, sz, nb, n_blocks, temp); // clique is stored as a collection of rectangles - n_blocks = (cl - 1) / nb_ + 1; + n_blocks = (cl - 1) / nb + 1; Int64 schur_size{}; for (Int j = 0; j < n_blocks; ++j) { - const Int jb = std::min(nb_, cl - j * nb_); - schur_size += (Int64)(cl - j * nb_) * jb; + const Int jb = std::min(nb, cl - j * nb); + schur_size += (Int64)(cl - j * nb) * jb; } cl_entries = schur_size; } @@ -1108,6 +1119,8 @@ void Analyse::reorderChildren() { } void Analyse::computeBlockStart() { + const Int nb = FH_opt_.nb; + clique_block_start_.resize(sn_count_); // compute starting position of each block of columns in the clique, for // each supernode @@ -1115,10 +1128,10 @@ void Analyse::computeBlockStart() { const Int sn_size = sn_start_[sn + 1] - sn_start_[sn]; const Int ldf = ptr_sn_[sn + 1] - ptr_sn_[sn]; const Int ldc = ldf - sn_size; - const Int n_blocks = (ldc - 1) / nb_ + 1; + const Int n_blocks = (ldc - 1) / nb + 1; Int64 schur_size = - getDiagStart(ldc, ldc, nb_, n_blocks, clique_block_start_[sn]); + getDiagStart(ldc, ldc, nb, n_blocks, clique_block_start_[sn]); clique_block_start_[sn].push_back(schur_size); } } @@ -1136,7 +1149,8 @@ Int Analyse::checkOverflow() const { const Int sn_size = sn_start_[sn + 1] - sn_start_[sn]; const Int front_size = ptr_sn_[sn + 1] - ptr_sn_[sn]; - if ((Int64)front_size * std::min(sn_size, nb_) > int32_limit) return 1; + if ((Int64)front_size * std::min(sn_size, FH_opt_.nb) > int32_limit) + return 1; } return 0; @@ -1261,7 +1275,6 @@ Int Analyse::run(Symbolic& S) { S.largest_front_ = *std::max_element(sn_indices_.begin(), sn_indices_.end()); S.serial_storage_ = serial_storage_; S.flops_ = dense_ops_; - S.block_size_ = nb_; S.max_stack_size_ = max_stack_size_; S.tree_depth_ = maxDepthTree(sn_parent_); diff --git a/highs/ipm/hipo/factorhighs/Analyse.h b/highs/ipm/hipo/factorhighs/Analyse.h index 330578615d6..a6432dafef6 100644 --- a/highs/ipm/hipo/factorhighs/Analyse.h +++ b/highs/ipm/hipo/factorhighs/Analyse.h @@ -5,6 +5,7 @@ #include #include "DataCollector.h" +#include "FactorHighsOptions.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" #include "ipm/hipo/auxiliary/Logger.h" @@ -79,8 +80,7 @@ class Analyse { std::vector> clique_block_start_{}; - // block size - Int nb_{}; + const FHoptions& FH_opt_; const Logger* logger_; DataCollector& data_; @@ -109,9 +109,9 @@ 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 Logger* logger, - DataCollector& data, const std::vector& perm); + Analyse(Int n, Int nz, const Int* rows, const Int* ptr, const Int* signs, + const FHoptions& FH_opt, const Logger* logger, DataCollector& data, + const Int* perm); // Run analyse phase and save the result in Symbolic object S Int run(Symbolic& S); diff --git a/highs/ipm/hipo/factorhighs/CallAndTimeBlas.cpp b/highs/ipm/hipo/factorhighs/CallAndTimeBlas.cpp index 6cabb4ff635..994178ab4de 100644 --- a/highs/ipm/hipo/factorhighs/CallAndTimeBlas.cpp +++ b/highs/ipm/hipo/factorhighs/CallAndTimeBlas.cpp @@ -2,7 +2,7 @@ #include "DataCollector.h" #include "DenseFact.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "HighsExternalApi.h" #include "Timing.h" #include "ipm/hipo/auxiliary/Auxiliary.h" diff --git a/highs/ipm/hipo/factorhighs/DataCollector.cpp b/highs/ipm/hipo/factorhighs/DataCollector.cpp index c0e70d48014..2f17bd13c14 100644 --- a/highs/ipm/hipo/factorhighs/DataCollector.cpp +++ b/highs/ipm/hipo/factorhighs/DataCollector.cpp @@ -1,6 +1,6 @@ #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "ipm/hipo/auxiliary/Logger.h" namespace hipo { diff --git a/highs/ipm/hipo/factorhighs/DataCollector.h b/highs/ipm/hipo/factorhighs/DataCollector.h index e5967787d8e..75593969e46 100644 --- a/highs/ipm/hipo/factorhighs/DataCollector.h +++ b/highs/ipm/hipo/factorhighs/DataCollector.h @@ -6,7 +6,7 @@ #include #include -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "Timing.h" #include "ipm/hipo/auxiliary/IntConfig.h" #include "ipm/hipo/auxiliary/Logger.h" diff --git a/highs/ipm/hipo/factorhighs/DenseFact.h b/highs/ipm/hipo/factorhighs/DenseFact.h index 5e12d09aac5..408a75d0d72 100644 --- a/highs/ipm/hipo/factorhighs/DenseFact.h +++ b/highs/ipm/hipo/factorhighs/DenseFact.h @@ -2,7 +2,8 @@ #define FACTORHIGHS_DENSE_FACT_H #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsOptions.h" +#include "FactorHighsSettings.h" #include "ipm/hipo/auxiliary/IntConfig.h" namespace hipo { @@ -31,14 +32,14 @@ namespace hipo { // dense factorisation kernel Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, - double thresh, const Regul& regul, double* totalreg, Int* swaps, - double* pivot_2x2, DataCollector& data); + double thresh, double* totalreg, Int* swaps, double* pivot_2x2, + DataCollector& data, const FHoptions& options); // dense partial factorisation, in "hybrid formats" -Int denseFactFH(char format, Int n, Int k, Int nb, double* A, double* B, - const Int* pivot_sign, double thresh, const Regul& regul, - double* totalreg, Int* swaps, double* pivot_2x2, bool parnode, - DataCollector& data); +Int denseFactFH(char format, Int n, Int k, double* A, double* B, + const Int* pivot_sign, double thresh, double* totalreg, + Int* swaps, double* pivot_2x2, bool parnode, + DataCollector& data, const FHoptions& options); // function to convert A from lower packed, to lower-blocked-hybrid format Int denseFactFP2FH(double* A, Int nrow, Int ncol, Int nb, DataCollector& data); diff --git a/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp b/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp index 666703f3cf3..ec0ba3e309e 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactHybrid.cpp @@ -10,10 +10,10 @@ namespace hipo { // Factorisation with "hybrid formats". -Int denseFactFH(char format, Int n, Int k, Int nb, double* A, double* B, - const Int* pivot_sign, double thresh, const Regul& regval, - double* totalreg, Int* swaps, double* pivot_2x2, bool parnode, - DataCollector& data) { +Int denseFactFH(char format, Int n, Int k, double* A, double* B, + const Int* pivot_sign, double thresh, double* totalreg, + Int* swaps, double* pivot_2x2, bool parnode, + DataCollector& data, const FHoptions& options) { // =========================================================================== // Partial blocked factorisation // Matrix A is in format FH @@ -30,6 +30,8 @@ Int denseFactFH(char format, Int n, Int k, Int nb, double* A, double* B, // quick return if (n == 0) return kRetOk; + const Int nb = options.nb; + // number of blocks of columns const Int n_blocks = (k - 1) / nb + 1; @@ -93,18 +95,18 @@ Int denseFactFH(char format, Int n, Int k, Int nb, double* A, double* B, &pivot_sign[j * nb] + jb); Int* swaps_current = &swaps[j * nb]; double* pivot_2x2_current = &pivot_2x2[j * nb]; - Int info = - denseFactK('U', jb, D, jb, pivot_sign_current.data(), thresh, regval, - regul_current, swaps_current, pivot_2x2_current, data); + Int info = denseFactK('U', jb, D, jb, pivot_sign_current.data(), thresh, + regul_current, swaps_current, pivot_2x2_current, data, + options); if (info != 0) return info; -#ifdef HIPO_PIVOTING - // swap columns in R - applySwaps(swaps_current, M, jb, R, data); + if (options.pivoting) { + // swap columns in R + applySwaps(swaps_current, M, jb, R, data); - // unswap regularisation, to keep it with original ordering - permuteWithSwaps(regul_current, swaps_current, jb, true); -#endif + // unswap regularisation, to keep it with original ordering + permuteWithSwaps(regul_current, swaps_current, jb, true); + } if (M > 0) { // =========================================================================== @@ -182,7 +184,7 @@ Int denseFactFH(char format, Int n, Int k, Int nb, double* A, double* B, offset += jb * col_jj; } HIPO_CLOCK_STOP(2, data, kTimeDenseFact_main); - + // =========================================================================== // UPDATE SCHUR COMPLEMENT // =========================================================================== diff --git a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp index eada1cdcf1d..cc3e0a9f34e 100644 --- a/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp +++ b/highs/ipm/hipo/factorhighs/DenseFactKernel.cpp @@ -1,7 +1,8 @@ #include "CallAndTimeBlas.h" #include "DataCollector.h" #include "DenseFact.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsOptions.h" +#include "FactorHighsSettings.h" #include "ReturnValues.h" #include "Swaps.h" #include "ipm/hipo/auxiliary/Auxiliary.h" @@ -38,21 +39,23 @@ static std::pair maxInCol(Int j, Int n, Int m, double* A, return {r, maxval}; } -static void staticReg(double& pivot, Int sign, const Regul& regval, +static void staticReg(double& pivot, Int sign, const FHoptions& options, double& totalreg) { // apply static regularisation double old_pivot = pivot; + if (sign == 0) sign = pivot > 0 ? 1 : -1; if (sign > 0) - pivot += regval.primal; + pivot += options.reg_p; else - pivot -= regval.dual; + pivot -= options.reg_d; totalreg = pivot - old_pivot; } -static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, - Int* sign, double thresh, const Regul& regval, - double* totalreg, DataCollector& data) { +static bool blockBunchKaufman(const Int j, Int n, double* A, Int lda, + Int* swaps, Int* sign, double thresh, + const FHoptions& options, double* totalreg, + DataCollector& data) { // Perform Bunch-Kaufman pivoting within a block of the supernode (see Schenk, // Gartner, ETNA 2006). // It works only for upper triangular A. @@ -82,24 +85,25 @@ static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, // Max in column j of diagonal block auto res = maxInCol(j, n, j, A, lda); double gamma_j = res.second; - Int r = res.first; - double Ajj = sign[j] > 0 ? A[j + lda * j] + regval.dual - : A[j + lda * j] - regval.primal; + const Int r = res.first; + const Int sign_j = sign[j] != 0 ? sign[j] : (A[j + lda * j] > 0 ? 1 : -1); + double Ajj = sign_j > 0 ? A[j + lda * j] + options.reg_d + : A[j + lda * j] - options.reg_p; - if (std::max(std::abs(Ajj), gamma_j) <= thresh || sign[j] * Ajj < 0 || + if (std::max(std::abs(Ajj), gamma_j) <= thresh || sign_j * Ajj < 0 || j == n - 1) { // Must accept current pivot double old_pivot = A[j + lda * j]; - staticReg(A[j + lda * j], sign[j], regval, totalreg[j]); + staticReg(A[j + lda * j], sign_j, options, totalreg[j]); - if (sign[j] * A[j + lda * j] < 0) { + if (sign_j * A[j + lda * j] < 0) { data.setWrongSign(A[j + lda * j]); - // A[j + lda * j] = sign[j] * thresh; + // A[j + lda * j] = sign_j * thresh; } if (std::max(std::abs(Ajj), gamma_j) < thresh) { // perturbe pivot - A[j + lda * j] = sign[j] * thresh; + A[j + lda * j] = sign_j * thresh; data.countRegPiv(); } totalreg[j] = A[j + lda * j] - old_pivot; @@ -110,34 +114,38 @@ static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, assert(r >= 0); res = maxInCol(j, n, r, A, lda); double gamma_r = res.second; - double Arr = sign[r] > 0 ? A[r + lda * r] + regval.primal - : A[r + lda * r] - regval.dual; + const Int sign_r = sign[r] != 0 ? sign[r] : (A[r + lda * r] > 0 ? 1 : -1); + double Arr = sign_r > 0 ? A[r + lda * r] + options.reg_p + : A[r + lda * r] - options.reg_d; if ((std::abs(Ajj) >= kAlphaBK * gamma_j || std::abs(Ajj) * gamma_r >= kAlphaBK * gamma_j * gamma_j)) { // Accept current pivot - staticReg(A[j + lda * j], sign[j], regval, totalreg[j]); + staticReg(A[j + lda * j], sign_j, options, totalreg[j]); - if (sign[j] * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); + if (sign_j * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); } else if (std::abs(Arr) >= kAlphaBK * gamma_r) { // Use pivot r swapCols('U', n, A, lda, j, r, swaps, sign, data); - staticReg(A[j + lda * j], sign[j], regval, totalreg[j]); + staticReg(A[j + lda * j], sign_j, options, totalreg[j]); - if (sign[j] * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); + if (sign_j * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); } else { // Use 2x2 pivot (j,r) + const Int sign_jp1 = sign[j + 1] != 0 + ? sign[j + 1] + : (A[j + 1 + lda * (j + 1)] > 0 ? 1 : -1); swapCols('U', n, A, lda, j + 1, r, swaps, sign, data); flag_2x2 = true; - staticReg(A[j + lda * j], sign[j], regval, totalreg[j]); - staticReg(A[j + 1 + lda * (j + 1)], sign[j + 1], regval, totalreg[j + 1]); + staticReg(A[j + lda * j], sign_j, options, totalreg[j]); + staticReg(A[j + 1 + lda * (j + 1)], sign_jp1, options, totalreg[j + 1]); - if (sign[j] * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); - if (sign[j + 1] * A[j + 1 + lda * (j + 1)] < 0) + if (sign_j * A[j + lda * j] < 0) data.setWrongSign(A[j + lda * j]); + if (sign_jp1 * A[j + 1 + lda * (j + 1)] < 0) data.setWrongSign(A[j + 1 + lda * (j + 1)]); } } @@ -147,8 +155,8 @@ static bool blockBunchKaufman(Int j, Int n, double* A, Int lda, Int* swaps, } Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, - double thresh, const Regul& regval, double* totalreg, Int* swaps, - double* pivot_2x2, DataCollector& data) { + double thresh, double* totalreg, Int* swaps, double* pivot_2x2, + DataCollector& data, const FHoptions& options) { // =========================================================================== // Factorisation kernel // Matrix A is in format F @@ -180,10 +188,10 @@ Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, for (Int j = 0; j < n; j += step) { bool flag_2x2 = false; -#ifdef HIPO_PIVOTING - flag_2x2 = blockBunchKaufman(j, n, A, lda, swaps, pivot_sign, thresh, - regval, totalreg, data); -#endif + if (options.pivoting) { + flag_2x2 = blockBunchKaufman(j, n, A, lda, swaps, pivot_sign, thresh, + options, totalreg, data); + } // cannot do 2x2 pivoting on last column assert(j < n - 1 || flag_2x2 == false); @@ -198,13 +206,13 @@ Int denseFactK(char uplo, Int n, double* A, Int lda, Int* pivot_sign, if (std::isnan(Ajj)) return kRetInvalidPivot; // add regularisation - staticReg(Ajj, pivot_sign[j], regval, totalreg[j]); + staticReg(Ajj, pivot_sign[j], options, totalreg[j]); -#ifndef HIPO_PIVOTING - // add static regularisation - staticReg(Ajj, pivot_sign[j], regul[j]); - data.setMaxReg(std::abs(regul[j])); -#endif + if (!options.pivoting) { + // add static regularisation + staticReg(Ajj, pivot_sign[j], options, totalreg[j]); + data.setMaxReg(std::abs(totalreg[j])); + } // save reciprocal of pivot A[j + lda * j] = 1.0 / Ajj; diff --git a/highs/ipm/hipo/factorhighs/DgemmParallel.cpp b/highs/ipm/hipo/factorhighs/DgemmParallel.cpp index 3ba954a6fdf..2cab7dcdcee 100644 --- a/highs/ipm/hipo/factorhighs/DgemmParallel.cpp +++ b/highs/ipm/hipo/factorhighs/DgemmParallel.cpp @@ -1,7 +1,7 @@ #include "DgemmParallel.h" #include "CallAndTimeBlas.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "parallel/HighsParallel.h" namespace hipo { diff --git a/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp b/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp deleted file mode 100644 index 3773532f8e6..00000000000 --- a/highs/ipm/hipo/factorhighs/FactorHiGHS.cpp +++ /dev/null @@ -1,57 +0,0 @@ -#include "FactorHiGHS.h" - -#include "Analyse.h" -#include "DataCollector.h" -#include "Factorise.h" -#include "ipm/hipo/auxiliary/Logger.h" - -namespace hipo { - -FHsolver::FHsolver(const Logger* logger, Int block_size) - : logger_{logger}, nb_{block_size > 0 ? block_size : default_nb_} { -#ifdef HIPO_COLLECT_EXPENSIVE_DATA - if (logger_) - logger_->printw( - "Running in debug mode: COLLECTING EXPENSIVE FACTORISATION DATA\n"); -#endif -#if HIPO_TIMING_LEVEL > 0 - if (logger_) - logger_->printw("Running in debug mode: COLLECTING EXPENSIVE TIMING DATA\n"); -#endif -} - -FHsolver::~FHsolver() { - if (logger_) { - data_.printTimes(*logger_); - data_.printIter(*logger_); - } -} - -void FHsolver::newIter() { data_.append(); } - -void FHsolver::setRegularisation(double reg_p, double reg_d) { - regul_.primal = reg_p; - regul_.dual = reg_d; -} - -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_, 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_, logger_, data_, sn_columns_, - &serial_stack_); - return fact_obj.run(N_); -} - -Int FHsolver::solve(std::vector& x) { return N_.solve(x); } - -void FHsolver::getRegularisation(std::vector& reg) { N_.getReg(reg); } - -} // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHiGHS.h b/highs/ipm/hipo/factorhighs/FactorHiGHS.h deleted file mode 100644 index 79a325db07a..00000000000 --- a/highs/ipm/hipo/factorhighs/FactorHiGHS.h +++ /dev/null @@ -1,110 +0,0 @@ -#ifndef FACTOR_HIGHS_H -#define FACTOR_HIGHS_H - -#include "CliqueStack.h" -#include "DataCollector.h" -#include "Numeric.h" -#include "Symbolic.h" -#include "ipm/hipo/auxiliary/Logger.h" - -/* - -Direct solver for IPM matrices. - -Consider a sparse symmetric matrix M in CSC format. -Only its lower triangular part is used; entries in the upper triangle are -ignored. -The matrix has n rows/cols and nz nonzero entries in the lower triangle. -It is stored using three arrays: -- ptr, column pointers, of length n+1; -- rows, row indices, of length nz; -- vals, values, of length nz. - -The direct solver uses the following objects: -- Symbolic, to store the symbolic factorization; -- FHsolver, to perform analyse and factorise phases. - -Define a vector signs that contains the expected sign of each pivot (1 or -1). -Define a right-hand side rhs, which will be overwritten with the solution of -M^{-1} * rhs. The pre-computed fill-reducing ordering to use is stored in the -vector perm. - -Then, the factorization is performed as follows. - - Symbolic S; - FHsolver FH; - FH.analyse(S, rows, ptr, signs, perm); - FH.factorise(S, rows, ptr, val); - FH.solve(x); - -Printing to screen is achieved using the interface in auxiliary/Logger.h. - ... - Logger logger; - FHsolver FH(&logger); - ... -Pass nothing to suppress all logging. - -To add static regularisation when the pivots are selected, use -setRegularisation(reg_p,reg_d) to choose values of primal and dual -regularisation. If regularisation is already added to the matrix, ignore. - -The default block size is 128. To set a different block size, pass it as second -input to the constructor. - -*/ - -namespace hipo { - -class FHsolver { - const Logger* logger_; - DataCollector data_; - Regul regul_; - Numeric N_; - CliqueStack serial_stack_; - - const Int nb_; // block size - static const Int default_nb_ = 128; - - // Columns of factorisation, stored by supernode. - // This memory is allocated the first time that it is used. Subsequent - // factorisations reuse the same memory. - std::vector> sn_columns_; - - public: - // Create object and initialise DataCollector - FHsolver(const Logger* logger = nullptr, Int block_size = default_nb_); - - // Print collected data (if any) and terminate DataCollector - ~FHsolver(); - - // Perform analyse phase of matrix with sparsity pattern given by rows and - // ptr, and store symbolic factorisation in object S. - // See ReturnValues.h for errors. - Int analyse(Symbolic& S, const std::vector& rows, - const std::vector& ptr, const std::vector& signs, - const std::vector& perm); - - // Perform factorise phase of matrix given by rows, ptr, vals, and store - // numerical factorisation in object N. See ReturnValues.h for errors. - Int factorise(const Symbolic& S, const std::vector& rows, - const std::vector& ptr, const std::vector& vals); - - // Perform solve phase with rhs given by x, which is overwritten with the - // solution. - Int solve(std::vector& x); - - // If multiple factorisation are performed, call newIter() before each - // factorisation. This is used only to collect data for debugging, if - // expensive data collection is turned on at compile time. - void newIter(); - - // Set values for static regularisation to be added when a pivot is selected. - // If regularisation is already added to the matrix, ignore. - void setRegularisation(double reg_p, double reg_d); - - void getRegularisation(std::vector& reg); -}; - -} // namespace hipo - -#endif \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHighs.cpp b/highs/ipm/hipo/factorhighs/FactorHighs.cpp new file mode 100644 index 00000000000..e50b4f24865 --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighs.cpp @@ -0,0 +1,195 @@ +#include "FactorHighs.h" + +#include "Analyse.h" +#include "DataCollector.h" +#include "Factorise.h" +#include "HighsExternalApi.h" +#include "ipm/hipo/auxiliary/Auxiliary.h" +#include "ipm/hipo/auxiliary/Logger.h" + +namespace hipo { + +FHsolver::~FHsolver() { + if (logger_) { + data_.printTimes(*logger_); + data_.printIter(*logger_); + } + if (local_logger_ && logger_) delete logger_; +} + +void FHsolver::newIter() { data_.append(); } + +void FHsolver::setRegularisation(double reg_p, double reg_d) { + options_.reg_p = reg_p; + options_.reg_d = reg_d; +} + +Int FHsolver::analyse(Symbolic& S, Int n, Int nz, const Int* rows, + const Int* ptr, const Int* signs, const Int* perm) { + Analyse an_obj(n, nz, rows, ptr, signs, options_, logger_, data_, perm); + return an_obj.run(S); +} + +Int FHsolver::factorise(const Symbolic& S, Int n, Int nz, const Int* rows, + const Int* ptr, const double* vals) { + Factorise fact_obj(S, n, nz, rows, ptr, vals, options_, logger_, data_, + sn_columns_, &serial_stack_); + return fact_obj.run(N_); +} + +Int FHsolver::solve(double* x, Int k) const { return N_.solve(x, k); } +Int FHsolver::forwardSolve(double* x, Int k) const { + return N_.forwardSolve(x, k); +} +Int FHsolver::diagSolve(double* x, Int k) const { return N_.diagSolve(x, k); } +Int FHsolver::backwardSolve(double* x, Int k) const { + return N_.backwardSolve(x, k); +} + +void FHsolver::getRegularisation(double* reg) { N_.getReg(reg); } + +void FHsolver::setBlockSize(Int nb) { + if (nb > 0) options_.nb = nb; +} + +void FHsolver::setPivoting(bool pivoting) { options_.pivoting = pivoting; } + +void FHsolver::setLogger(const Logger* logger, bool use_printf) { + if (local_logger_ && logger_) delete logger_; + local_logger_ = false; + + logger_ = logger; + if (!logger_) { + logger_ = new Logger(use_printf); + if (logger_) local_logger_ = true; + } + +#ifdef HIPO_COLLECT_EXPENSIVE_DATA + if (logger_) + logger_->printw( + "Running in debug mode: COLLECTING EXPENSIVE FACTORISATION DATA\n"); +#endif +#if HIPO_TIMING_LEVEL > 0 + if (logger_) + logger_->printw( + "Running in debug mode: COLLECTING EXPENSIVE TIMING DATA\n"); +#endif +} + +void FHsolver::inertia(Int& pos, Int& neg, Int& zero, double tol) const { + N_.inertia(pos, neg, zero, tol); +} + +void FHsolver::setOneIndexing(bool one_indexing) { + options_.one_indexing = one_indexing; +} + +void FHsolver::iperm(const Symbolic& S, Int* ip) const { + const std::vector& iperm = S.iperm(); + const Int offset = options_.one_indexing ? 1 : 0; + for (Int i = 0; i < static_cast(iperm.size()); ++i) { + ip[i] = iperm[i] + offset; + } +} + +static void getFull(Int n, Int nz, const Int* rows, const Int* ptr, + bool one_indexing, std::vector& full_ptr_v, + std::vector& full_rows_v) { + std::vector rows_tmp(nz); + std::vector ptr_tmp(n + 1); + for (Int i = 0; i < nz; ++i) rows_tmp[i] = rows[i]; + for (Int i = 0; i < n + 1; ++i) ptr_tmp[i] = ptr[i]; + + if (one_indexing) { + for (Int& i : rows_tmp) --i; + for (Int& i : ptr_tmp) --i; + } + + // metis, amd, rcm require a full copy of the matrix, without diagonal + // entries + fullFromLower(n, nz, ptr_tmp.data(), rows_tmp.data(), full_ptr_v, + full_rows_v); +} + +Int FHsolver::reorderMetis(Int n, Int nz, const Int* rows, const Int* ptr, + Int* perm, bool full_matrix_0) const { + const Int *full_ptr, *full_rows; + std::vector full_ptr_v, full_rows_v; + if (full_matrix_0) { + full_ptr = ptr; + full_rows = rows; + } else { + getFull(n, nz, rows, ptr, options_.one_indexing, full_ptr_v, full_rows_v); + full_ptr = full_ptr_v.data(); + full_rows = full_rows_v.data(); + } + + idx_t options[METIS_NOPTIONS]; + HighsExtras::metis::set_default_options(options); + options[METIS_OPTION_SEED] = kMetisSeed; + + options[METIS_OPTION_DBGLVL] = 0; + + // no2hop improves the quality of ordering in general + options[METIS_OPTION_NO2HOP] = 1; + + std::vector iperm(n); + + Int status = HighsExtras::metis::nodeND(&n, full_ptr, full_rows, nullptr, + options, perm, iperm.data()); + + if (options_.one_indexing) + for (Int i = 0; i < n; ++i) perm[i]++; + + return status != METIS_OK; +} + +Int FHsolver::reorderAmd(Int n, Int nz, const Int* rows, const Int* ptr, + Int* perm, bool full_matrix_0) const { + const Int *full_ptr, *full_rows; + std::vector full_ptr_v, full_rows_v; + if (full_matrix_0) { + full_ptr = ptr; + full_rows = rows; + } else { + getFull(n, nz, rows, ptr, options_.one_indexing, full_ptr_v, full_rows_v); + full_ptr = full_ptr_v.data(); + full_rows = full_rows_v.data(); + } + + double control[AMD_CONTROL]; + HighsExtras::amd::set_defaults(control); + double info[AMD_INFO]; + + Int status = + HighsExtras::amd::order(n, full_ptr, full_rows, perm, control, info); + + if (options_.one_indexing) + for (Int i = 0; i < n; ++i) perm[i]++; + + return status != AMD_OK; +} + +Int FHsolver::reorderRcm(Int n, Int nz, const Int* rows, const Int* ptr, + Int* perm, bool full_matrix_0) const { + const Int *full_ptr, *full_rows; + std::vector full_ptr_v, full_rows_v; + if (full_matrix_0) { + full_ptr = ptr; + full_rows = rows; + } else { + getFull(n, nz, rows, ptr, options_.one_indexing, full_ptr_v, full_rows_v); + full_ptr = full_ptr_v.data(); + full_rows = full_rows_v.data(); + } + + Int status = + HighsExtras::rcm::genrcm(n, full_ptr[n], full_ptr, full_rows, perm); + + if (options_.one_indexing) + for (Int i = 0; i < n; ++i) perm[i]++; + + return status != 0; +} + +} // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHighs.h b/highs/ipm/hipo/factorhighs/FactorHighs.h new file mode 100644 index 00000000000..c83918bc3a0 --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighs.h @@ -0,0 +1,175 @@ +#ifndef FACTOR_HIGHS_H +#define FACTOR_HIGHS_H + +#include "CliqueStack.h" +#include "DataCollector.h" +#include "FactorHighsOptions.h" +#include "Numeric.h" +#include "Symbolic.h" +#include "ipm/hipo/auxiliary/Logger.h" + +/* + +Direct solver for IPM matrices. + +Consider a sparse symmetric matrix M in CSC format. +Only its lower triangular part is used; entries in the upper triangle are +ignored. Diagonal entries must be present for each column, even if the entry is +zero. +The matrix has n rows/cols and nz nonzero entries. It is stored using +three arrays: +- ptr, column pointers, of length n+1; +- rows, row indices, of length nz; +- vals, values, of length nz. + +The direct solver uses the following objects: +- Symbolic, to store the symbolic factorization; +- FHsolver, to perform analyse and factorise phases. + +Define a right-hand side rhs, which will be overwritten with the solution of +M^{-1} * rhs. Use reorderMetis, reorderAmd, reorderRcm to obtain a fill-reducing +ordering of the matrix, or use a precomputed one. The ordering is stored in the +array perm. + +Zero-based indexing is assumed for the sparsity pattern and permutation. Use +setOneIndexing(true) to use one-based indexing. + +Define a vector signs that contains the expected sign of each pivot: +- 1 for pivots expected to be positive +- -1 for pivots expected to be negative +- 0 for pivots without an expected sign. +This is used to determine the sign of the regularisation to apply. + +Then, the factorization is performed as follows. + + Symbolic S; + FHsolver FH; + FH.analyse(S, n, nz rows, ptr, signs, perm); + FH.factorise(S, n, nz rows, ptr, val); + FH.solve(x); + +Printing to screen is achieved using the interface in auxiliary/Logger.h. +Use setLogger to pass the Logger object to use: FH.setLogger(&logger). +To use printf instead, use FH.setLogger(nullptr,true). +Logging is off by default. Use FH.setLogger(nullptr, false) to switch logging +off. + +To add static regularisation when the pivots are selected, use +setRegularisation(reg_p,reg_d) to choose values of primal and dual +regularisation. If regularisation is already added to the matrix, ignore. + +The expected sign of a pivot is used to determine which static regularisation to +apply. Pivots with negative expected sign are regularised with -reg_p, pivots +with positive expected sign are regularised with +reg_d. If the sign is unknown, +the computed sign of the pivot is used to determine which regularisation to +apply. + +Notice that this sign convention corresponds to a saddle point system with +negative (1,1) block and positive (2,2) block. + +Dynamic regularisation may perturb a pivot if it is too small, and it +does so based on the expected sign of the pivot. This procedure uses the +computed sign of the pivot, if the sign is unknown. + +Notice that the fill-reducing ordering can be modified during the call to +analyse. The inverse permutation used during the factorisation can be accessed +via the Symbolic object, S.iperm(). + +*/ + +namespace hipo { + +class FHsolver { + const Logger* logger_; + DataCollector data_; + Numeric N_; + CliqueStack serial_stack_; + bool local_logger_ = false; + + FHoptions options_{}; + + // Columns of factorisation, stored by supernode. + // This memory is allocated the first time that it is used. Subsequent + // factorisations reuse the same memory. + std::vector> sn_columns_; + + public: + // Print collected data (if any) and terminate DataCollector + ~FHsolver(); + + // Perform analyse phase of matrix with sparsity pattern given by rows and + // ptr, and store symbolic factorisation in object S. + // See ReturnValues.h for errors. + Int analyse(Symbolic& S, Int n, Int nz, const Int* rows, const Int* ptr, + const Int* signs, const Int* perm); + + // Perform factorise phase of matrix given by rows, ptr, vals, and store + // numerical factorisation in object N. See ReturnValues.h for errors. + Int factorise(const Symbolic& S, Int n, Int nz, const Int* rows, + const Int* ptr, const double* vals); + + // Perform solve phase with rhs given by x, which is overwritten with the + // solution. Multiple rhs are supported by doing the solves in parallel. + Int solve(double* x, Int k = 1) const; + + // Perform partial solves with L, D, L^T (including permutation) + Int forwardSolve(double* x, Int k = 1) const; + Int diagSolve(double* x, Int k = 1) const; + Int backwardSolve(double* x, Int k = 1) const; + + // If multiple factorisation are performed, call newIter() before each + // factorisation. This is used only to collect data for debugging, if + // expensive data collection is turned on at compile time. + void newIter(); + + // Set values for static regularisation to be added when a pivot is selected. + // If regularisation is already added to the matrix, ignore. + void setRegularisation(double reg_p = 0.0, double reg_d = 0.0); + + // Extract the regularisation values used, including static and dynamic. + void getRegularisation(double* reg); + + // Set block size for dense linear algebra + void setBlockSize(Int nb); + + // Set pivoting optins, on by default. + // It uses a static variation of Bunch-Kaufman pivoting, with potential + // dynamic regularisation. If pivoting is switched off, only static + // regularisation is applied. + void setPivoting(bool pivoting); + + // Pass the Logger object to be used for logging. Alternatively, printf can be + // used for logging, by passing a nullptr and setting use_printf to true. + // By default, logging is off. + void setLogger(const Logger* logger = nullptr, bool use_printf = false); + const Logger* getLogger() const { return logger_; } + + // Compute number of positive, negative and zero pivots, using tol as + // tolerance for zero. + void inertia(Int& pos, Int& neg, Int& zero, double tol = 1e-16) const; + + // Set options for 1-based indexing, off by default. + // If set to true, the vectors passed to analyse and factorise are assumed to + // use 1-based indexing, so that all entries are shifted down by 1. + void setOneIndexing(bool one_indexing); + + void iperm(const Symbolic& S, Int* ip) const; + + // Compute fill-reducing permutation using nested dissection (Metis), approx + // minimum degree (amd) or reverse Cuthill-McGee (rcm). + // Set full_matrix_0 = false if the sparsity pattern in rows and ptr + // corresponds to the lower triangle of the matrix. + // Set full_matrix_0 = true if the sparsity pattern in rows and ptr + // corresponds to the full matrix (without the diagonal entries), with 0-based + // indexing. + Int reorderMetis(Int n, Int nz, const Int* rows, const Int* ptr, Int* perm, + bool full_matrix_0) const; + Int reorderAmd(Int n, Int nz, const Int* rows, const Int* ptr, Int* perm, + bool full_matrix_0) const; + Int reorderRcm(Int n, Int nz, const Int* rows, const Int* ptr, Int* perm, + bool full_matrix_0) const; +}; + +} // namespace hipo + +#endif \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHighsOptions.h b/highs/ipm/hipo/factorhighs/FactorHighsOptions.h new file mode 100644 index 00000000000..a937e43d244 --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighsOptions.h @@ -0,0 +1,19 @@ +#ifndef FACTOR_HIGHS_OPTIONS_H +#define FACTOR_HIGHS_OPTIONS_H + +#include "FactorHighsSettings.h" +#include "ipm/hipo/auxiliary/IntConfig.h" + +namespace hipo { + +struct FHoptions { + double reg_p = 0.0; + double reg_d = 0.0; + Int nb = kBlockSize; + bool pivoting = true; + bool one_indexing = false; +}; + +} // namespace hipo + +#endif \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h b/highs/ipm/hipo/factorhighs/FactorHighsSettings.h similarity index 87% rename from highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h rename to highs/ipm/hipo/factorhighs/FactorHighsSettings.h index 81725710add..e3cb9983a4b 100644 --- a/highs/ipm/hipo/factorhighs/FactorHiGHSSettings.h +++ b/highs/ipm/hipo/factorhighs/FactorHighsSettings.h @@ -9,11 +9,6 @@ // SWITCHES // =========================================================================== -// Switch on/off pivoting. It uses a static variation of Bunch-Kaufman pivoting, -// with potential dynamic regularisation. If pivoting is switched off, only -// static regularisation is applied. -#define HIPO_PIVOTING - // Collect data during regularisation, e.g. number of regularised pivots, 2x2 // pivots, pivot swaps, pivots with wrong sign, min and max entry of L and D. // This can be quite expensive and should only be used for debugging. @@ -32,6 +27,8 @@ namespace hipo { +const Int kBlockSize = 128; + // supernode amalgamation const Int kStartThreshRelax = 256; const double kUpperRatioRelax = 0.02; diff --git a/highs/ipm/hipo/factorhighs/FactorHighs_c_api.cpp b/highs/ipm/hipo/factorhighs/FactorHighs_c_api.cpp new file mode 100644 index 00000000000..80fbbec900d --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighs_c_api.cpp @@ -0,0 +1,105 @@ +#include "FactorHighs_c_api.h" + +#include "FactorHighs.h" +#include "HighsExternalApi.h" +#include "parallel/HighsParallel.h" + +HighsInt FactorHighs_initialise(HighsInt threads) { + highs::parallel::initialize_scheduler(threads); + bool hipo_available = HighsExternalApi::isAvailable(); + if (hipo_available) HighsExtras::blas::openblas_set_num_threads(1); + return !hipo_available; +} + +void FactorHighs_terminate(void) { HighsTaskExecutor::shutdown(true); } + +void* FactorHighs_create(void) { return new hipo::FHsolver(); } +void FactorHighs_destroy(void* FH) { delete (hipo::FHsolver*)FH; } + +void* FactorHighs_symbolic_create(void) { return new hipo::Symbolic(); } +void FactorHighs_symbolic_destroy(void* S) { delete (hipo::Symbolic*)S; } + +HighsInt FactorHighs_analyse(void* FH, void* S, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + const HighsInt* signs, const HighsInt* perm) { + return ((hipo::FHsolver*)FH) + ->analyse(*(hipo::Symbolic*)S, n, nz, rows, ptr, signs, perm); +} + +HighsInt FactorHighs_factorise(void* FH, const void* S, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + const double* vals) { + return ((hipo::FHsolver*)FH) + ->factorise(*(hipo::Symbolic*)S, n, nz, rows, ptr, vals); +} + +HighsInt FactorHighs_solve(void* FH, double* x, HighsInt k) { + return ((hipo::FHsolver*)FH)->solve(x, k); +} + +HighsInt FactorHighs_forwardSolve(void* FH, double* x, HighsInt k) { + return ((hipo::FHsolver*)FH)->forwardSolve(x, k); +} + +HighsInt FactorHighs_diagSolve(void* FH, double* x, HighsInt k) { + return ((hipo::FHsolver*)FH)->diagSolve(x, k); +} + +HighsInt FactorHighs_backwardSolve(void* FH, double* x, HighsInt k) { + return ((hipo::FHsolver*)FH)->backwardSolve(x, k); +} + +void FactorHighs_setRegularisation(void* FH, double reg_p, double reg_d) { + ((hipo::FHsolver*)FH)->setRegularisation(reg_p, reg_d); +} + +void FactorHighs_getRegularisation(void* FH, double* reg) { + ((hipo::FHsolver*)FH)->getRegularisation(reg); +} + +void FactorHighs_newIter(void* FH) { ((hipo::FHsolver*)FH)->newIter(); } + +void FactorHighs_setBlockSize(void* FH, HighsInt nb) { + ((hipo::FHsolver*)FH)->setBlockSize(nb); +} + +void FactorHighs_setPivoting(void* FH, HighsInt pivoting) { + ((hipo::FHsolver*)FH)->setPivoting(pivoting); +} + +void FactorHighs_setLogging(void* FH, HighsInt display) { + ((hipo::FHsolver*)FH)->setLogger(nullptr, display); +} + +void FactorHighs_inertia(void* FH, HighsInt* pos, HighsInt* neg, HighsInt* zero, + double tol) { + ((hipo::FHsolver*)FH)->inertia(*pos, *neg, *zero, tol); +} + +void FactorHighs_symbolic_print(void* FH, void* S, HighsInt verbose) { + ((hipo::Symbolic*)S)->print(*(((hipo::FHsolver*)FH)->getLogger()), verbose); +} + +void FactorHighs_iperm(void* FH, void* S, HighsInt* ip) { + ((hipo::FHsolver*)FH)->iperm(*(hipo::Symbolic*)S, ip); +} + +void FactorHighs_setOneIndexing(void* FH, HighsInt one_indexing) { + ((hipo::FHsolver*)FH)->setOneIndexing(one_indexing); +} + +HighsInt FactorHighs_reorderMetis(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm) { + return ((hipo::FHsolver*)FH)->reorderMetis(n, nz, rows, ptr, perm, 0); +} +HighsInt FactorHighs_reorderAmd(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm) { + return ((hipo::FHsolver*)FH)->reorderAmd(n, nz, rows, ptr, perm, 0); +} +HighsInt FactorHighs_reorderRcm(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm) { + return ((hipo::FHsolver*)FH)->reorderRcm(n, nz, rows, ptr, perm, 0); +} \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHighs_c_api.h b/highs/ipm/hipo/factorhighs/FactorHighs_c_api.h new file mode 100644 index 00000000000..10c060b8de0 --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighs_c_api.h @@ -0,0 +1,70 @@ +#ifndef FACTOR_HIGHS_C_API_H +#define FACTOR_HIGHS_C_API_H + +#include "util/HighsInt.h" + +/* C API to HiPO linear solver + It is meant to be used outside of HiGHS as a standalone linear solver. + Refer to FactorHighs.h for a description of the functions. + Refer to FactorHighs_c_api_example.c for a small example. +*/ + +#ifdef __cplusplus +extern "C" { +#endif + +// Initialise and create/destroy objects +HighsInt FactorHighs_initialise(HighsInt threads); +void FactorHighs_terminate(void); +void* FactorHighs_symbolic_create(void); +void FactorHighs_symbolic_destroy(void* S); +void* FactorHighs_create(void); +void FactorHighs_destroy(void* FH); + +// Main functions +HighsInt FactorHighs_analyse(void* FH, void* S, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + const HighsInt* signs, const HighsInt* perm); +HighsInt FactorHighs_factorise(void* FH, const void* S, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + const double* vals); +HighsInt FactorHighs_solve(void* FH, double* x, HighsInt k); + +// Partial solves +HighsInt FactorHighs_forwardSolve(void* FH, double* x, HighsInt k); +HighsInt FactorHighs_diagSolve(void* FH, double* x, HighsInt k); +HighsInt FactorHighs_backwardSolve(void* FH, double* x, HighsInt k); + +// Reordering +HighsInt FactorHighs_reorderMetis(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm); +HighsInt FactorHighs_reorderAmd(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm); +HighsInt FactorHighs_reorderRcm(void* FH, HighsInt n, HighsInt nz, + const HighsInt* rows, const HighsInt* ptr, + HighsInt* perm); + +// Set options +void FactorHighs_setRegularisation(void* FH, double reg_p, double reg_d); +void FactorHighs_setBlockSize(void* FH, HighsInt nb); +void FactorHighs_setPivoting(void* FH, HighsInt pivoting); +void FactorHighs_setLogging(void* FH, HighsInt display); +void FactorHighs_setOneIndexing(void* FH, HighsInt one_indexing); + +// Extract information +void FactorHighs_getRegularisation(void* FH, double* reg); +void FactorHighs_iperm(void* FH, void* S, HighsInt* ip); +void FactorHighs_inertia(void* FH, HighsInt* pos, HighsInt* neg, HighsInt* zero, + double tol); + +// Other +void FactorHighs_symbolic_print(void* FH, void* S, HighsInt verbose); +void FactorHighs_newIter(void* FH); + +#ifdef __cplusplus +} +#endif + +#endif \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/FactorHighs_c_api_example.c b/highs/ipm/hipo/factorhighs/FactorHighs_c_api_example.c new file mode 100644 index 00000000000..3ae96509f10 --- /dev/null +++ b/highs/ipm/hipo/factorhighs/FactorHighs_c_api_example.c @@ -0,0 +1,143 @@ +#include "ipm/hipo/factorhighs/FactorHighs_c_api.h" +#include "math.h" +#include "stdio.h" + +/* + This example factorises the following symmetric matrix + + 5 + 0 3 + 3 2 9 + 4 0 -1 8 + 0 0 1 0 1 + + and solves a linear system, using the FactorHighs C API. + This file needs to be linked with HiGHS: + + clang FactorHighs_c_api_example.c + -I/path/to/HiGHS/highs -I/path/to/HiGHS/build + -L/path/to/HiGHS/build/lib -lhighs + -Wl,-rpath,/path/to/HiGHS/build/lib + +*/ + +int main() { + // problem size + enum { N = 5, NZ = 10, NP1 = N + 1 }; + + // define the lower triangle in CSC format, with 0-based indexing + int ptr[NP1] = {0, 3, 5, 8, 9, 10}; + int rows[NZ] = {0, 2, 3, 1, 2, 2, 3, 4, 3, 4}; + double vals[NZ] = {5, 3, 4, 3, 2, 9, -1, 1, 8, 1}; + + /* + 1-based indexing is available as well. + We define the problem shifted by 1: + + int ptr[n + 1] = {1, 4, 6, 9, 10, 11}; + int rows[nz] = {1, 3, 4, 2, 3, 3, 4, 5, 4, 5}; + double vals[nz] = {5, 3, 4, 3, 2, 9, -1, 1, 8, 1}; + + and set the corresponding option: + + FactorHighs_setOneIndexing(FH, 1); + */ + + // matrix is spd, so expect all pivots to be positive + int signs[N] = {1, 1, 1, 1, 1}; + + // rhs and expected solution + double rhs[N] = {1, 2, 3, 4, 5}; + const double lhs[N] = {0.457627118644068, 1.118644067796610, + -0.677966101694915, 0.186440677966102, + 5.677966101694915}; + + // initialise with default number of threads + const int num_threads = 0; + int initialise_status = FactorHighs_initialise(num_threads); + if (initialise_status) { + printf("HiPO or its dependencies are not available in this build\n"); + return 1; + } + void* S = FactorHighs_symbolic_create(); + void* FH = FactorHighs_create(); + + // set options + const int logging_on = 1; + FactorHighs_setLogging(FH, logging_on); + + const int block_size = 64; + FactorHighs_setBlockSize(FH, block_size); + + const int pivoting_off = 0; + FactorHighs_setPivoting(FH, pivoting_off); + + FactorHighs_setRegularisation(FH, 0.0, 0.0); + + // compute ordering with metis + int perm[N]; + int metis_status = FactorHighs_reorderMetis(FH, N, NZ, rows, ptr, perm); + if (metis_status) { + printf("Failed to compute a fill-reducing ordering with Metis\n"); + return 1; + } + + // perform analyse phase + int analyse_status = + FactorHighs_analyse(FH, S, N, NZ, rows, ptr, signs, perm); + if (analyse_status) { + printf("Failed to compute symbolic factorisation\n"); + return 1; + } + + // print extended statistics of symbolic factorisation + int verbose = 1; + FactorHighs_symbolic_print(FH, S, verbose); + + // print inverse permutation, that may have been modified by analyse phase + int iperm[N]; + FactorHighs_iperm(FH, S, iperm); + printf("\niperm: "); + for (int i = 0; i < N; ++i) printf("%d ", iperm[i]); + printf("\n"); + + // factorise the matrix + int factorise_status = FactorHighs_factorise(FH, S, N, NZ, rows, ptr, vals); + if (factorise_status) { + printf("Failed to factorise the matrix\n"); + return 1; + } + + // compute the inertia of the factorisation + int pos, neg, zero; + double zero_tolerance = 1e-16; + FactorHighs_inertia(FH, &pos, &neg, &zero, zero_tolerance); + printf("\nInertia (+/-/0): %d %d %d\n", pos, neg, zero); + + // triangular solve + int solve_status = FactorHighs_solve(FH, rhs, 1); + if (solve_status) { + printf("Failed to solve the linear system\n"); + return 1; + } + + /* + could also perform: + + FactorHighs_forwardSolve(FH, rhs, 1); + FactorHighs_diagSolve(FH, rhs, 1); + FactorHighs_backwardSolve(FH, rhs, 1); + */ + + // compute error + double error = 0.0; + for (int i = 0; i < N; ++i) error += fabs(rhs[i] - lhs[i]); + printf("\nError: %e\n", error); + + // terminate + FactorHighs_symbolic_destroy(S); + FactorHighs_destroy(FH); + FactorHighs_terminate(); + + return 0; +} diff --git a/highs/ipm/hipo/factorhighs/Factorise.cpp b/highs/ipm/hipo/factorhighs/Factorise.cpp index 0ffeb5ec71b..bb2d40c0ac0 100644 --- a/highs/ipm/hipo/factorhighs/Factorise.cpp +++ b/highs/ipm/hipo/factorhighs/Factorise.cpp @@ -4,7 +4,7 @@ #include #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "FormatHandler.h" #include "HybridHybridFormatHandler.h" #include "ReturnValues.h" @@ -14,23 +14,24 @@ namespace hipo { -Factorise::Factorise(const Symbolic& S, const std::vector& rowsM, - const std::vector& ptrM, - const std::vector& valM, const Regul& regul, - const Logger* logger, DataCollector& data, +Factorise::Factorise(const Symbolic& S, Int n, Int nz, const Int* rowsM, + const Int* ptrM, const double* valM, + const FHoptions& FH_opt, const Logger* logger, + DataCollector& data, std::vector>& sn_columns, CliqueStack* stack) : S_{S}, sn_columns_{sn_columns}, - regul_{regul}, logger_{logger}, data_{data}, + FH_opt_{FH_opt}, stack_{stack} { // Input the symmetric matrix to be factorised in CSC format and the symbolic // factorisation coming from Analyse. // Only the lower triangular part of the matrix is used. + // Diagonal entries must be stored for each column, even if zero. - n_ = ptrM.size() - 1; + n_ = n; if (n_ != S_.size()) { if (logger_) @@ -41,9 +42,15 @@ Factorise::Factorise(const Symbolic& S, const std::vector& rowsM, } // Make a copy of the matrix to be factorised - rowsM_ = rowsM; - valM_ = valM; - ptrM_ = ptrM; + rowsM_ = std::vector(rowsM, rowsM + nz); + valM_ = std::vector(valM, valM + nz); + ptrM_ = std::vector(ptrM, ptrM + n + 1); + + // adjust data if one-based indexing is used + if (FH_opt_.one_indexing) { + for (Int& i : rowsM_) --i; + for (Int& i : ptrM_) --i; + } // Permute the matrix. // This also removes any entry not in the lower triangle. @@ -141,7 +148,7 @@ void Factorise::processSupernode(Int sn) { // initialise the format handler // this also allocates space for the frontal matrix and schur complement std::unique_ptr FH(new HybridHybridFormatHandler( - S_, sn, regul_, data_, sn_columns_[sn], clique_ptr)); + S_, sn, data_, sn_columns_[sn], clique_ptr, FH_opt_)); HIPO_CLOCK_STOP(2, data_, kTimeFactorisePrepare); @@ -185,7 +192,8 @@ void Factorise::processSupernode(Int sn) { child_clique = schur_contribution_[child_sn].data(); if (!child_clique) { - if (logger_) logger_->printInfo("Missing child supernode contribution\n"); + if (logger_) + logger_->printInfo("Missing child supernode contribution\n"); flag_stop_.store(true, std::memory_order_relaxed); return; } @@ -326,6 +334,8 @@ bool Factorise::run(Numeric& num) { if (flag_stop_.load(std::memory_order_relaxed)) return true; + permuteVector(total_reg_, S_.iperm()); + // move factorisation to numerical object num.S_ = &S_; num.sn_columns_ = &sn_columns_; @@ -333,6 +343,7 @@ bool Factorise::run(Numeric& num) { num.swaps_ = std::move(swaps_); num.pivot_2x2_ = std::move(pivot_2x2_); num.data_ = &data_; + num.options_ = &FH_opt_; HIPO_CLOCK_STOP(1, data_, kTimeFactorise); diff --git a/highs/ipm/hipo/factorhighs/Factorise.h b/highs/ipm/hipo/factorhighs/Factorise.h index 0a673b5c50f..46f3b27daf3 100644 --- a/highs/ipm/hipo/factorhighs/Factorise.h +++ b/highs/ipm/hipo/factorhighs/Factorise.h @@ -4,6 +4,7 @@ #include #include "CliqueStack.h" +#include "FactorHighsOptions.h" #include "Numeric.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" @@ -56,14 +57,12 @@ class Factorise { // regularisation std::vector total_reg_{}; - // values for static regularisation - const Regul& regul_; - // flag to stop computation std::atomic flag_stop_{false}; const Logger* logger_; DataCollector& data_; + const FHoptions& FH_opt_; CliqueStack* stack_; @@ -72,10 +71,10 @@ class Factorise { void processSupernode(Int sn); public: - Factorise(const Symbolic& S, const std::vector& rowsM, - const std::vector& ptrM, const std::vector& valM, - const Regul& regul, const Logger* logger, DataCollector& data, - std::vector>& sn_columns, CliqueStack* stack); + Factorise(const Symbolic& S, Int n, Int nz, const Int* rowsM, const Int* ptrM, + const double* valM, const FHoptions& FH_opt, const Logger* logger, + DataCollector& data, std::vector>& sn_columns, + CliqueStack* stack); bool run(Numeric& num); }; diff --git a/highs/ipm/hipo/factorhighs/FormatHandler.cpp b/highs/ipm/hipo/factorhighs/FormatHandler.cpp index c3171fdd61b..ee07e12fcfc 100644 --- a/highs/ipm/hipo/factorhighs/FormatHandler.cpp +++ b/highs/ipm/hipo/factorhighs/FormatHandler.cpp @@ -6,15 +6,15 @@ namespace hipo { -FormatHandler::FormatHandler(const Symbolic& S, Int sn, const Regul& regul, - std::vector& frontal, double* clique_ptr) +FormatHandler::FormatHandler(const Symbolic& S, Int sn, + std::vector& frontal, double* clique_ptr, + const FHoptions& FH_opt) : S_{&S}, - regul_{regul}, sn_{sn}, - nb_{S_->blockSize()}, sn_size_{S_->snStart(sn_ + 1) - S_->snStart(sn_)}, ldf_{(Int)(S_->ptr(sn_ + 1) - S_->ptr(sn_))}, ldc_{ldf_ - sn_size_}, + FH_opt_{FH_opt}, frontal_{frontal}, clique_ptr_{clique_ptr} { local_reg_.resize(sn_size_); diff --git a/highs/ipm/hipo/factorhighs/FormatHandler.h b/highs/ipm/hipo/factorhighs/FormatHandler.h index 54f18a0a504..5b26415951e 100644 --- a/highs/ipm/hipo/factorhighs/FormatHandler.h +++ b/highs/ipm/hipo/factorhighs/FormatHandler.h @@ -3,7 +3,8 @@ #include -#include "FactorHiGHSSettings.h" +#include "FactorHighsOptions.h" +#include "FactorHighsSettings.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" @@ -30,14 +31,9 @@ class FormatHandler { // symbolic object const Symbolic* S_; - const Regul& regul_; - // supernode being processed const Int sn_{}; - // block size - const Int nb_{}; - // size of the supernode const Int sn_size_{}; @@ -47,6 +43,8 @@ class FormatHandler { // size of the clique const Int ldc_{}; + const FHoptions& FH_opt_; + // local copies to be moved at the end std::vector& frontal_; std::vector clique_{}; @@ -56,8 +54,8 @@ class FormatHandler { std::vector pivot_2x2_{}; public: - FormatHandler(const Symbolic& S, Int sn, const Regul& regul, - std::vector& frontal, double* clique_ptr); + FormatHandler(const Symbolic& S, Int sn, std::vector& frontal, + double* clique_ptr, const FHoptions& FH_opt); void terminate(std::vector& clique, std::vector& total_reg, std::vector& swaps, std::vector& pivot_2x2); diff --git a/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp b/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp index 4b17afdb77c..6144fb58afa 100644 --- a/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp +++ b/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp @@ -11,9 +11,11 @@ namespace hipo { HybridHybridFormatHandler::HybridHybridFormatHandler( - const Symbolic& S, Int sn, const Regul& regul, DataCollector& data, - std::vector& frontal, double* clique_ptr) - : FormatHandler(S, sn, regul, frontal, clique_ptr), data_{data} { + const Symbolic& S, Int sn, DataCollector& data, + std::vector& frontal, double* clique_ptr, const FHoptions& FH_opt) + : FormatHandler(S, sn, frontal, clique_ptr, FH_opt), + data_{data}, + nb_{FH_opt_.nb} { // initialise frontal and clique initFrontal(); @@ -94,10 +96,9 @@ Int HybridHybridFormatHandler::denseFactorise(double reg_thresh) { Int sn_start = S_->snStart(sn_); const Int* pivot_sign = &S_->pivotSign().data()[sn_start]; - status = denseFactFH('H', ldf_, sn_size_, S_->blockSize(), frontal_.data(), - clique_ptr_, pivot_sign, reg_thresh, regul_, - local_reg_.data(), swaps_.data(), pivot_2x2_.data(), - S_->parNode(), data_); + status = denseFactFH('H', ldf_, sn_size_, frontal_.data(), clique_ptr_, + pivot_sign, reg_thresh, local_reg_.data(), swaps_.data(), + pivot_2x2_.data(), S_->parNode(), data_, FH_opt_); return status; } diff --git a/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.h b/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.h index dea4772381c..ea811942340 100644 --- a/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.h +++ b/highs/ipm/hipo/factorhighs/HybridHybridFormatHandler.h @@ -9,6 +9,7 @@ namespace hipo { class HybridHybridFormatHandler : public FormatHandler { std::vector diag_start_; DataCollector& data_; + const Int nb_; void initFrontal() override; void initClique() override; @@ -21,9 +22,9 @@ class HybridHybridFormatHandler : public FormatHandler { void extremeEntries() override; public: - HybridHybridFormatHandler(const Symbolic& S, Int sn, const Regul& regul, - DataCollector& data, std::vector& frontal, - double* clique_ptr); + HybridHybridFormatHandler(const Symbolic& S, Int sn, DataCollector& data, + std::vector& frontal, double* clique_ptr, + const FHoptions& FH_opt); }; } // namespace hipo diff --git a/highs/ipm/hipo/factorhighs/HybridSolveHandler.cpp b/highs/ipm/hipo/factorhighs/HybridSolveHandler.cpp index 0fc6c5a6834..7d8603aa138 100644 --- a/highs/ipm/hipo/factorhighs/HybridSolveHandler.cpp +++ b/highs/ipm/hipo/factorhighs/HybridSolveHandler.cpp @@ -2,7 +2,7 @@ #include "CallAndTimeBlas.h" #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "FormatHandler.h" #include "Swaps.h" #include "ipm/hipo/auxiliary/Auxiliary.h" @@ -12,10 +12,13 @@ namespace hipo { HybridSolveHandler::HybridSolveHandler( const Symbolic& S, const std::vector>& sn_columns, const std::vector>& swaps, - const std::vector>& pivot_2x2, DataCollector& data) - : SolveHandler(S, sn_columns, data), swaps_{swaps}, pivot_2x2_{pivot_2x2} {} + const std::vector>& pivot_2x2, DataCollector& data, + const FHoptions& options) + : SolveHandler(S, sn_columns, data, options), + swaps_{swaps}, + pivot_2x2_{pivot_2x2} {} -void HybridSolveHandler::forwardSolve(std::vector& x) const { +void HybridSolveHandler::forwardSolve(double* x) const { // Forward solve. // Blas calls: dtrsv, dgemv @@ -23,7 +26,7 @@ void HybridSolveHandler::forwardSolve(std::vector& x) const { HIPO_CLOCK_CREATE; - const Int nb = S_.blockSize(); + const Int nb = options_.nb; for (Int sn = 0; sn < S_.sn(); ++sn) { // leading size of supernode @@ -50,13 +53,13 @@ void HybridSolveHandler::forwardSolve(std::vector& x) const { const Int jb = sn_size; const Int x_start = sn_start; -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply swaps to portion of rhs that is affected const Int* current_swaps = swaps_[sn].data(); - permuteWithSwaps(&x[x_start], current_swaps, jb); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply swaps to portion of rhs that is affected + permuteWithSwaps(&x[x_start], current_swaps, jb); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } HIPO_CLOCK_START(2); for (Int row = 0; row < jb; ++row) { @@ -74,12 +77,12 @@ void HybridSolveHandler::forwardSolve(std::vector& x) const { } HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_dense); -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply inverse swaps - permuteWithSwaps(&x[x_start], current_swaps, jb, true); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply inverse swaps + permuteWithSwaps(&x[x_start], current_swaps, jb, true); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } } else { // go through blocks of columns for this supernode @@ -93,13 +96,13 @@ void HybridSolveHandler::forwardSolve(std::vector& x) const { // index to access vector x const Int x_start = sn_start + nb * j; -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply swaps to portion of rhs that is affected const Int* current_swaps = &swaps_[sn][nb * j]; - permuteWithSwaps(&x[x_start], current_swaps, jb); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply swaps to portion of rhs that is affected + permuteWithSwaps(&x[x_start], current_swaps, jb); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } HIPO_CLOCK_START(2); callAndTime_dtrsv('U', 'T', 'U', jb, &sn_columns_[sn][SnCol_ind], jb, @@ -127,18 +130,18 @@ void HybridSolveHandler::forwardSolve(std::vector& x) const { HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_sparse); } -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply inverse swaps - permuteWithSwaps(&x[x_start], current_swaps, jb, true); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply inverse swaps + permuteWithSwaps(&x[x_start], current_swaps, jb, true); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } } } } } -void HybridSolveHandler::backwardSolve(std::vector& x) const { +void HybridSolveHandler::backwardSolve(double* x) const { // Backward solve. // Blas calls: dtrsv, dgemv @@ -146,7 +149,7 @@ void HybridSolveHandler::backwardSolve(std::vector& x) const { HIPO_CLOCK_CREATE; - const Int nb = S_.blockSize(); + const Int nb = options_.nb; // go through the sn in reverse order for (Int sn = S_.sn() - 1; sn >= 0; --sn) { @@ -175,13 +178,13 @@ void HybridSolveHandler::backwardSolve(std::vector& x) const { const Int jb = sn_size; const Int x_start = sn_start; -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply swaps to portion of rhs that is affected const Int* current_swaps = swaps_[sn].data(); - permuteWithSwaps(&x[x_start], current_swaps, jb); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply swaps to portion of rhs that is affected + permuteWithSwaps(&x[x_start], current_swaps, jb); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } HIPO_CLOCK_START(2); for (Int row = ldSn - 1; row >= jb; --row) { @@ -199,12 +202,12 @@ void HybridSolveHandler::backwardSolve(std::vector& x) const { } HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_dense); -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply inverse swaps - permuteWithSwaps(&x[x_start], current_swaps, jb, true); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply inverse swaps + permuteWithSwaps(&x[x_start], current_swaps, jb, true); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } } else { // go through blocks of columns for this supernode in reverse order @@ -218,13 +221,13 @@ void HybridSolveHandler::backwardSolve(std::vector& x) const { // index to access vector x const Int x_start = sn_start + nb * j; -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply swaps to portion of rhs that is affected const Int* current_swaps = &swaps_[sn][nb * j]; - permuteWithSwaps(&x[x_start], current_swaps, jb); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply swaps to portion of rhs that is affected + permuteWithSwaps(&x[x_start], current_swaps, jb); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } // temporary space for gemv const Int gemv_space = ldSn - nb * j - jb; @@ -252,25 +255,25 @@ void HybridSolveHandler::backwardSolve(std::vector& x) const { &x[x_start], 1, data_); HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_dense); -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply inverse swaps - permuteWithSwaps(&x[x_start], current_swaps, jb, true); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply inverse swaps + permuteWithSwaps(&x[x_start], current_swaps, jb, true); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } } } } } -void HybridSolveHandler::diagSolve(std::vector& x) const { +void HybridSolveHandler::diagSolve(double* x) const { // Diagonal solve // supernode columns in format FH HIPO_CLOCK_CREATE; - const Int nb = S_.blockSize(); + const Int nb = options_.nb; for (Int sn = 0; sn < S_.sn(); ++sn) { // leading size of supernode @@ -293,13 +296,13 @@ void HybridSolveHandler::diagSolve(std::vector& x) const { // number of columns in the block const Int jb = std::min(nb, sn_size - nb * j); -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply swaps to portion of rhs that is affected const Int* current_swaps = &swaps_[sn][nb * j]; - permuteWithSwaps(&x[sn_start + nb * j], current_swaps, jb); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply swaps to portion of rhs that is affected + permuteWithSwaps(&x[sn_start + nb * j], current_swaps, jb); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } HIPO_CLOCK_START(2); @@ -333,12 +336,107 @@ void HybridSolveHandler::diagSolve(std::vector& x) const { HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_dense); -#ifdef HIPO_PIVOTING - HIPO_CLOCK_START(2); - // apply inverse swaps - permuteWithSwaps(&x[sn_start + nb * j], current_swaps, jb, true); - HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); -#endif + if (options_.pivoting) { + HIPO_CLOCK_START(2); + // apply inverse swaps + permuteWithSwaps(&x[sn_start + nb * j], current_swaps, jb, true); + HIPO_CLOCK_STOP(2, data_, kTimeSolveSolve_swap); + } + + // move diag_start forward by number of diagonal entries in block + diag_start += jb * jb; + + // move diag_start forward by number of sub-diagonal entries in block + diag_start += (ldSn - nb * j - jb) * jb; + } + } +} + +void HybridSolveHandler::inertia(Int& pos, Int& neg, Int& zero, + double tol) const { + pos = 0; + neg = 0; + zero = 0; + + const Int nb = options_.nb; + + for (Int sn = 0; sn < S_.sn(); ++sn) { + // leading size of supernode + const Int ldSn = S_.ptr(sn + 1) - S_.ptr(sn); + + // number of columns in the supernode + const Int sn_size = S_.snStart(sn + 1) - S_.snStart(sn); + + // number of blocks of columns + const Int n_blocks = (sn_size - 1) / nb + 1; + + // index to access diagonal part of block + Int diag_start{}; + + // go through blocks of columns for this supernode + for (Int j = 0; j < n_blocks; ++j) { + // number of columns in the block + const Int jb = std::min(nb, sn_size - nb * j); + const double* current_2x2 = &pivot_2x2_[sn][nb * j]; + Int step = 1; + + // go through columns of block + for (Int col = 0; col < jb; col += step) { + if (current_2x2[col] == 0.0) { + // 1x1 pivots + step = 1; + const double inv_d = sn_columns_[sn][diag_start + col + jb * col]; + const double pivot = 1.0 / inv_d; + + if (pivot > tol) + ++pos; + else if (pivot < -tol) + ++neg; + else + ++zero; + + } else { + // 2x2 pivots + step = 2; + + // inverse of 2x2 pivot + const double i_d1 = sn_columns_[sn][diag_start + col + jb * col]; + const double i_d2 = + sn_columns_[sn][diag_start + col + 1 + jb * (col + 1)]; + const double i_off = current_2x2[col]; + + // determinant and trace of inverse + const double i_det = i_d1 * i_d2 - i_off * i_off; + const double i_trace = i_d1 + i_d2; + + // determinant and trace of 2x2 pivot + const double det = 1.0 / i_det; + const double trace = det * i_trace; + + if (std::abs(det) < tol) { + // det is zero, so at least one pivot is zero + ++zero; + + if (trace > tol) + ++pos; + else if (trace < -tol) + ++neg; + else + ++zero; + + } else if (det > tol) { + // det is positive, so pivots have same sign + if (trace > 0) + pos += 2; + else + neg -= 2; + } else { + // indefinite 2x2 pivot + ++pos; + ++neg; + } + } + } // move diag_start forward by number of diagonal entries in block diag_start += jb * jb; diff --git a/highs/ipm/hipo/factorhighs/HybridSolveHandler.h b/highs/ipm/hipo/factorhighs/HybridSolveHandler.h index e490946fbaf..5ed2d7871d0 100644 --- a/highs/ipm/hipo/factorhighs/HybridSolveHandler.h +++ b/highs/ipm/hipo/factorhighs/HybridSolveHandler.h @@ -10,15 +10,16 @@ class HybridSolveHandler : public SolveHandler { const std::vector>& pivot_2x2_; public: - void forwardSolve(std::vector& x) const override; - void backwardSolve(std::vector& x) const override; - void diagSolve(std::vector& x) const override; + void forwardSolve(double* x) const override; + void backwardSolve(double* x) const override; + void diagSolve(double* x) const override; + void inertia(Int& pos, Int& neg, Int& zero, double tol) const override; HybridSolveHandler(const Symbolic& S, const std::vector>& sn_columns, const std::vector>& swaps, const std::vector>& pivot_2x2, - DataCollector& data); + DataCollector& data, const FHoptions& options); }; } // namespace hipo diff --git a/highs/ipm/hipo/factorhighs/KrylovMethodsIpm.cpp b/highs/ipm/hipo/factorhighs/KrylovMethodsIpm.cpp index 1a59826a626..7a7c57acda5 100644 --- a/highs/ipm/hipo/factorhighs/KrylovMethodsIpm.cpp +++ b/highs/ipm/hipo/factorhighs/KrylovMethodsIpm.cpp @@ -54,7 +54,7 @@ void IpmMatrix::apply(std::vector& x) const { } IpmFactor::IpmFactor(const Numeric& N) : N_{N} {} -void IpmFactor::apply(std::vector& x) const { N_.solve(x); } +void IpmFactor::apply(std::vector& x) const { N_.solve(x.data()); } void NeDiagPrec::reset(const HighsSparseMatrix& A, const std::vector& scaling) { diff --git a/highs/ipm/hipo/factorhighs/Numeric.cpp b/highs/ipm/hipo/factorhighs/Numeric.cpp index 5872d5d9d3f..0330bfe1c39 100644 --- a/highs/ipm/hipo/factorhighs/Numeric.cpp +++ b/highs/ipm/hipo/factorhighs/Numeric.cpp @@ -1,7 +1,7 @@ #include "Numeric.h" #include "DataCollector.h" -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" #include "HybridSolveHandler.h" #include "ReturnValues.h" #include "Timing.h" @@ -12,15 +12,16 @@ namespace hipo { -Int Numeric::solve(std::vector& x) const { +Int Numeric::solve(double* x) const { // Return the number of solves performed - if (!sn_columns_ || !S_) return kRetInvalidPointer; + if (!sn_columns_ || !S_ || !data_ || !options_) return kRetInvalidPointer; HIPO_CLOCK_CREATE; // initialise solve handler - HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_); + HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_, + *options_); // permute rhs HIPO_CLOCK_START(2); @@ -43,11 +44,65 @@ Int Numeric::solve(std::vector& x) const { return kRetOk; } -void Numeric::getReg(std::vector& reg) { - // unpermute regularisation - permuteVector(total_reg_, S_->iperm()); +Int Numeric::forwardSolve(double* x) const { + if (!sn_columns_ || !S_ || !data_ || !options_) return kRetInvalidPointer; + HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_, + *options_); + permuteVectorInverse(x, S_->iperm()); + SH.forwardSolve(x); + return kRetOk; +} +Int Numeric::diagSolve(double* x) const { + if (!sn_columns_ || !S_ || !data_ || !options_) return kRetInvalidPointer; + HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_, + *options_); + SH.diagSolve(x); + return kRetOk; +} +Int Numeric::backwardSolve(double* x) const { + if (!sn_columns_ || !S_ || !data_ || !options_) return kRetInvalidPointer; + HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_, + *options_); + SH.backwardSolve(x); + permuteVector(x, S_->iperm()); + return kRetOk; +} + +#define SOLVE_MULTIPLE(f) \ + if (k == 1) \ + return f(x); \ + else { \ + highs::parallel::TaskGroup tg; \ + const Int n = S_->size(); \ + std::atomic fail{false}; \ + for (Int i = 0; i < k; ++i) { \ + tg.spawn([=, &fail]() { \ + Int status = f(&x[i * n]); \ + if (status) fail.store(true, std::memory_order_relaxed); \ + }); \ + } \ + tg.taskWait(); \ + return fail; \ + } + +Int Numeric::solve(double* x, Int k) const { SOLVE_MULTIPLE(solve); } +Int Numeric::forwardSolve(double* x, Int k) const { + SOLVE_MULTIPLE(forwardSolve); +} +Int Numeric::diagSolve(double* x, Int k) const { SOLVE_MULTIPLE(diagSolve); } +Int Numeric::backwardSolve(double* x, Int k) const { + SOLVE_MULTIPLE(backwardSolve); +} + +void Numeric::getReg(double* reg) { + std::memcpy(reg, total_reg_.data(), total_reg_.size() * sizeof(double)); +} - reg = std::move(total_reg_); +void Numeric::inertia(Int& pos, Int& neg, Int& zero, double tol) const { + if (!sn_columns_ || !S_ || !data_ || !options_) return; + HybridSolveHandler SH(*S_, *sn_columns_, swaps_, pivot_2x2_, *data_, + *options_); + SH.inertia(pos, neg, zero, tol); } } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/Numeric.h b/highs/ipm/hipo/factorhighs/Numeric.h index 0480f050126..a647b201a78 100644 --- a/highs/ipm/hipo/factorhighs/Numeric.h +++ b/highs/ipm/hipo/factorhighs/Numeric.h @@ -5,6 +5,7 @@ #include #include "DataCollector.h" +#include "FactorHighsOptions.h" #include "SolveHandler.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" @@ -30,15 +31,27 @@ class Numeric { DataCollector* data_ = nullptr; + const FHoptions* options_; + friend class Factorise; // dynamic regularisation applied to the matrix std::vector total_reg_{}; public: - // Full solve - Int solve(std::vector& x) const; - void getReg(std::vector& reg); + Int solve(double* x) const; + Int solve(double* x, Int k) const; + + Int forwardSolve(double* x) const; + Int diagSolve(double* x) const; + Int backwardSolve(double* x) const; + + Int forwardSolve(double* x, Int k) const; + Int diagSolve(double* x, Int k) const; + Int backwardSolve(double* x, Int k) const; + + void getReg(double* reg); + void inertia(Int& pos, Int& neg, Int& zero, double tol) const; }; } // namespace hipo diff --git a/highs/ipm/hipo/factorhighs/SolveHandler.cpp b/highs/ipm/hipo/factorhighs/SolveHandler.cpp index 7ad20cc95f9..8c95f16f34c 100644 --- a/highs/ipm/hipo/factorhighs/SolveHandler.cpp +++ b/highs/ipm/hipo/factorhighs/SolveHandler.cpp @@ -4,7 +4,7 @@ namespace hipo { SolveHandler::SolveHandler(const Symbolic& S, const std::vector>& sn_columns, - DataCollector& data) - : S_{S}, sn_columns_{sn_columns}, data_{data} {} + DataCollector& data, const FHoptions& options) + : S_{S}, sn_columns_{sn_columns}, data_{data}, options_{options} {} } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/factorhighs/SolveHandler.h b/highs/ipm/hipo/factorhighs/SolveHandler.h index 491f923be4f..524453a1a4c 100644 --- a/highs/ipm/hipo/factorhighs/SolveHandler.h +++ b/highs/ipm/hipo/factorhighs/SolveHandler.h @@ -6,6 +6,7 @@ #include "DataCollector.h" #include "Symbolic.h" #include "ipm/hipo/auxiliary/IntConfig.h" +#include "FactorHighsOptions.h" namespace hipo { @@ -22,10 +23,12 @@ class SolveHandler { const std::vector>& sn_columns_; DataCollector& data_; + const FHoptions& options_; + public: SolveHandler(const Symbolic& S, const std::vector>& sn_columns, - DataCollector& data); + DataCollector& data, const FHoptions& options); // avoid copies SolveHandler(const SolveHandler&) = delete; @@ -38,9 +41,11 @@ class SolveHandler { // Pure virtual functions. // These need to be defined by any derived class. // ================================================================= - virtual void forwardSolve(std::vector& x) const = 0; - virtual void backwardSolve(std::vector& x) const = 0; - virtual void diagSolve(std::vector& x) const = 0; + virtual void forwardSolve(double* x) const = 0; + virtual void backwardSolve(double* x) const = 0; + virtual void diagSolve(double* x) const = 0; + + virtual void inertia(Int& pos, Int& neg, Int& zero, double tol) const = 0; }; } // namespace hipo diff --git a/highs/ipm/hipo/factorhighs/Symbolic.cpp b/highs/ipm/hipo/factorhighs/Symbolic.cpp index 163154281d2..d5802c7a08e 100644 --- a/highs/ipm/hipo/factorhighs/Symbolic.cpp +++ b/highs/ipm/hipo/factorhighs/Symbolic.cpp @@ -2,7 +2,7 @@ #include -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" namespace hipo { @@ -17,7 +17,6 @@ Int64 Symbolic::nz() const { return nz_; } double Symbolic::flops() const { return flops_; } double Symbolic::spops() const { return spops_; } double Symbolic::critops() const { return critops_; } -Int Symbolic::blockSize() const { return block_size_; } Int Symbolic::size() const { return n_; } Int Symbolic::sn() const { return sn_; } double Symbolic::fillin() const { return fillin_; } diff --git a/highs/ipm/hipo/factorhighs/Symbolic.h b/highs/ipm/hipo/factorhighs/Symbolic.h index 1c1dff63d65..c0b9fbe575b 100644 --- a/highs/ipm/hipo/factorhighs/Symbolic.h +++ b/highs/ipm/hipo/factorhighs/Symbolic.h @@ -14,9 +14,6 @@ class Symbolic { bool parallel_tree_ = false; bool parallel_node_ = false; - // Size of blocks for dense factorisation - Int block_size_; - // Statistics about symbolic factorisation Int n_{}; Int64 nz_{}; @@ -109,7 +106,6 @@ class Symbolic { double flops() const; double spops() const; double critops() const; - Int blockSize() const; Int size() const; Int sn() const; double fillin() const; diff --git a/highs/ipm/hipo/factorhighs/Timing.h b/highs/ipm/hipo/factorhighs/Timing.h index 0a0a221a6fe..630aeae1df7 100644 --- a/highs/ipm/hipo/factorhighs/Timing.h +++ b/highs/ipm/hipo/factorhighs/Timing.h @@ -1,7 +1,7 @@ #ifndef FACTORHIGHS_TIMING_H #define FACTORHIGHS_TIMING_H -#include "FactorHiGHSSettings.h" +#include "FactorHighsSettings.h" namespace hipo { diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHighsSolver.cpp similarity index 88% rename from highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp rename to highs/ipm/hipo/ipm/FactorHighsSolver.cpp index 50a66ef846f..8c0624262b3 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp +++ b/highs/ipm/hipo/ipm/FactorHighsSolver.cpp @@ -1,8 +1,7 @@ -#include "FactorHiGHSSolver.h" +#include "FactorHighsSolver.h" #include -#include "HighsExternalApi.h" #include "Status.h" #include "ipm/hipo/auxiliary/Auxiliary.h" #include "ipm/hipo/auxiliary/Logger.h" @@ -10,11 +9,11 @@ namespace hipo { -FactorHiGHSSolver::FactorHiGHSSolver(KktMatrix& kkt, Options& options, +FactorHighsSolver::FactorHighsSolver(KktMatrix& kkt, Options& options, const Model& model, const Regularisation& regul, Info& info, IpmData& record, const Logger& logger) - : FH_(&logger, options.block_size), + : FH_{}, S_{}, kkt_{kkt}, regul_{regul}, @@ -22,9 +21,12 @@ FactorHiGHSSolver::FactorHiGHSSolver(KktMatrix& kkt, Options& options, data_{record}, logger_{logger}, model_{model}, - options_{options} {} + options_{options} { + FH_.setBlockSize(options.block_size); + FH_.setLogger(&logger); +} -void FactorHiGHSSolver::clear() { +void FactorHighsSolver::clear() { valid_ = false; FH_.newIter(); } @@ -33,7 +35,7 @@ void FactorHiGHSSolver::clear() { // Analyse phase // ========================================================================= -Int FactorHiGHSSolver::analyseAS(Symbolic& S) { +Int FactorHighsSolver::analyseAS(Symbolic& S) { // Perform analyse phase of augmented system and return symbolic factorisation // in object S and the status. @@ -56,7 +58,7 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) { return status; } -Int FactorHiGHSSolver::analyseNE(Symbolic& S) { +Int FactorHighsSolver::analyseNE(Symbolic& S) { // Perform analyse phase of normal equations and return symbolic factorisation // in object S and the status. Structure of the matrix must be already // computed. @@ -80,7 +82,7 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S) { // Factorise phase // ========================================================================= -Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { +Int FactorHighsSolver::factorAS(const std::vector& scaling) { // only execute factorisation if it has not been done yet assert(!this->valid_); @@ -90,7 +92,8 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { FH_.setRegularisation(regul_.primal, regul_.dual); Clock clock; - if (FH_.factorise(S_, kkt_.rowsAS, kkt_.ptrAS, kkt_.valAS)) + if (FH_.factorise(S_, kkt_.n(), kkt_.nz(), kkt_.rowsAS.data(), + kkt_.ptrAS.data(), kkt_.valAS.data())) return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -99,7 +102,7 @@ Int FactorHiGHSSolver::factorAS(const std::vector& scaling) { return kStatusOk; } -Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { +Int FactorHighsSolver::factorNE(const std::vector& scaling) { // only execute factorisation if it has not been done yet assert(!this->valid_); @@ -109,7 +112,8 @@ Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { FH_.setRegularisation(regul_.primal, regul_.dual); Clock clock; - if (FH_.factorise(S_, kkt_.rowsNE, kkt_.ptrNE, kkt_.valNE)) + if (FH_.factorise(S_, kkt_.n(), kkt_.nz(), kkt_.rowsNE.data(), + kkt_.ptrNE.data(), kkt_.valNE.data())) return kStatusErrorFactorise; info_.factor_time += clock.stop(); info_.factor_number++; @@ -122,7 +126,7 @@ Int FactorHiGHSSolver::factorNE(const std::vector& scaling) { // Solve phase // ========================================================================= -Int FactorHiGHSSolver::solveAS(const std::vector& rhs_x, +Int FactorHighsSolver::solveAS(const std::vector& rhs_x, const std::vector& rhs_y, std::vector& lhs_x, std::vector& lhs_y) { @@ -137,7 +141,7 @@ Int FactorHiGHSSolver::solveAS(const std::vector& rhs_x, rhs.insert(rhs.end(), rhs_y.begin(), rhs_y.end()); Clock clock; - if (FH_.solve(rhs)) return kStatusErrorSolve; + if (FH_.solve(rhs.data())) return kStatusErrorSolve; info_.solve_time += clock.stop(); info_.solve_number++; @@ -151,7 +155,7 @@ Int FactorHiGHSSolver::solveAS(const std::vector& rhs_x, return kStatusOk; } -Int FactorHiGHSSolver::solveNE(const std::vector& rhs, +Int FactorHighsSolver::solveNE(const std::vector& rhs, std::vector& lhs) { // only execute the solve if factorisation is valid assert(this->valid_); @@ -160,7 +164,7 @@ Int FactorHiGHSSolver::solveNE(const std::vector& rhs, lhs = rhs; Clock clock; - if (FH_.solve(lhs)) return kStatusErrorSolve; + if (FH_.solve(lhs.data())) return kStatusErrorSolve; info_.solve_time += clock.stop(); info_.solve_number++; @@ -174,7 +178,7 @@ Int FactorHiGHSSolver::solveNE(const std::vector& rhs, // Automatic selections // ========================================================================= -Int FactorHiGHSSolver::setup() { +Int FactorHighsSolver::setup() { Clock clock; if (Int status = setNla()) return status; @@ -198,7 +202,7 @@ Int FactorHiGHSSolver::setup() { return kStatusOk; } -Int FactorHiGHSSolver::chooseNla() { +Int FactorHighsSolver::chooseNla() { // Choose whether to use augmented system or normal equations. Symbolic symb_NE{}; @@ -353,7 +357,7 @@ Int FactorHiGHSSolver::chooseNla() { return status; } -Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, +Int FactorHighsSolver::chooseOrdering(const std::vector& rows, const std::vector& ptr, const std::vector& signs, Symbolic& S, std::string& ordering, @@ -385,7 +389,8 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, // compute full-format matrix without diagonal entries std::vector full_ptr, full_rows; - fullFromLower(ptr, rows, full_ptr, full_rows); + fullFromLower(ptr.size() - 1, rows.size(), ptr.data(), rows.data(), full_ptr, + full_rows); Int n = full_ptr.size() - 1; std::vector perm(n), iperm(n); @@ -399,36 +404,22 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, nla.c_str()); if (orderings_to_try[i] == kHipoMetisString) { - idx_t options[METIS_NOPTIONS]; - HighsExtras::metis::set_default_options(options); - options[METIS_OPTION_SEED] = kMetisSeed; - - options[METIS_OPTION_DBGLVL] = 0; - - // no2hop improves the quality of ordering in general - options[METIS_OPTION_NO2HOP] = 1; - - std::vector iperm(n); - Int status = HighsExtras::metis::nodeND( - &n, full_ptr.data(), full_rows.data(), nullptr, options, - permutations[i].data(), iperm.data()); - if (status != METIS_OK) failure[i] = true; + Int status = + FH_.reorderMetis(n, rows.size(), full_rows.data(), full_ptr.data(), + permutations[i].data(), true); + if (status) failure[i] = true; } else if (orderings_to_try[i] == kHipoAmdString) { - double control[AMD_CONTROL]; - HighsExtras::amd::set_defaults(control); - double info[AMD_INFO]; - Int status = - HighsExtras::amd::order(n, full_ptr.data(), full_rows.data(), - permutations[i].data(), control, info); - if (status != AMD_OK) failure[i] = true; + FH_.reorderAmd(n, rows.size(), full_rows.data(), full_ptr.data(), + permutations[i].data(), true); + if (status) failure[i] = true; } else if (orderings_to_try[i] == kHipoRcmString) { Int status = - HighsExtras::rcm::genrcm(n, full_ptr.back(), full_ptr.data(), - full_rows.data(), permutations[i].data()); - if (status != 0) failure[i] = true; + FH_.reorderRcm(n, rows.size(), full_rows.data(), full_ptr.data(), + permutations[i].data(), true); + if (status) failure[i] = true; } else { assert(1 == 0); @@ -443,7 +434,9 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, return; } - failure[i] = FH_.analyse(symbolics[i], rows, ptr, signs, permutations[i]); + failure[i] = + FH_.analyse(symbolics[i], ptr.size() - 1, rows.size(), rows.data(), + ptr.data(), signs.data(), permutations[i].data()); if (failure[i] && logger_.debug(2)) { logger_.print("Failed symbolic:"); @@ -522,7 +515,7 @@ Int FactorHiGHSSolver::chooseOrdering(const std::vector& rows, return num_success > 0 ? kStatusOk : kStatusErrorAnalyse; } -Int FactorHiGHSSolver::setNla() { +Int FactorHighsSolver::setNla() { std::stringstream log_stream; if (options_.nla == kHipoNormalEqString && model_.nonSeparableQp()) { @@ -572,7 +565,7 @@ Int FactorHiGHSSolver::setNla() { return kStatusOk; } -void FactorHiGHSSolver::setParallel() { +void FactorHighsSolver::setParallel() { // Set parallel options bool parallel_tree = false; bool parallel_node = false; @@ -657,11 +650,11 @@ void FactorHiGHSSolver::setParallel() { // Other stuff // ========================================================================= -double FactorHiGHSSolver::flops() const { return S_.flops(); } -double FactorHiGHSSolver::spops() const { return S_.spops(); } -double FactorHiGHSSolver::nz() const { return (double)S_.nz(); } -void FactorHiGHSSolver::getReg(std::vector& reg) { - FH_.getRegularisation(reg); +double FactorHighsSolver::flops() const { return S_.flops(); } +double FactorHighsSolver::spops() const { return S_.spops(); } +double FactorHighsSolver::nz() const { return (double)S_.nz(); } +void FactorHighsSolver::getReg(std::vector& reg) { + FH_.getRegularisation(reg.data()); } } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHighsSolver.h similarity index 91% rename from highs/ipm/hipo/ipm/FactorHiGHSSolver.h rename to highs/ipm/hipo/ipm/FactorHighsSolver.h index 651147aadb2..c23b63f7f4f 100644 --- a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h +++ b/highs/ipm/hipo/ipm/FactorHighsSolver.h @@ -10,11 +10,11 @@ #include "LinearSolver.h" #include "Model.h" #include "ipm/hipo/auxiliary/IntConfig.h" -#include "ipm/hipo/factorhighs/FactorHiGHS.h" +#include "ipm/hipo/factorhighs/FactorHighs.h" namespace hipo { -class FactorHiGHSSolver : public LinearSolver { +class FactorHighsSolver : public LinearSolver { // object to perform factorisation FHsolver FH_; @@ -43,7 +43,7 @@ class FactorHiGHSSolver : public LinearSolver { Int analyseNE(Symbolic& S); public: - FactorHiGHSSolver(KktMatrix& kkt, Options& options, const Model& model, + FactorHighsSolver(KktMatrix& kkt, Options& options, const Model& model, const Regularisation& regul, Info& info, IpmData& record, const Logger& logger); diff --git a/highs/ipm/hipo/ipm/Iterate.cpp b/highs/ipm/hipo/ipm/Iterate.cpp index 84cbd355791..68a2663e337 100644 --- a/highs/ipm/hipo/ipm/Iterate.cpp +++ b/highs/ipm/hipo/ipm/Iterate.cpp @@ -2,7 +2,7 @@ #include "Parameters.h" #include "ipm/IpxWrapper.h" -#include "ipm/hipo/factorhighs/FactorHiGHSSettings.h" +#include "ipm/hipo/factorhighs/FactorHighsSettings.h" #include "model/HighsHessianUtils.h" namespace hipo { @@ -486,6 +486,11 @@ Int Iterate::finalResiduals(Info& info) const { } void Iterate::getReg(LinearSolver& LS, const std::string& nla) { + if (nla == kHipoNormalEqString) + total_reg.resize(model.m()); + else + total_reg.resize(model.m() + model.n()); + // extract regularisation LS.getReg(total_reg); diff --git a/highs/ipm/hipo/ipm/KktMatrix.cpp b/highs/ipm/hipo/ipm/KktMatrix.cpp index 304ea874ff0..840fdaa09d1 100644 --- a/highs/ipm/hipo/ipm/KktMatrix.cpp +++ b/highs/ipm/hipo/ipm/KktMatrix.cpp @@ -262,4 +262,20 @@ void KktMatrix::freeNEmemory() { freeVector(corr_A); } +Int KktMatrix::n() const { + if (nla() == kHipoNormalEqString) return ptrNE.size() - 1; + if (nla() == kHipoAugmentedString) return ptrAS.size() - 1; + return -1; +} +Int KktMatrix::nz() const { + if (nla() == kHipoNormalEqString) return rowsNE.size(); + if (nla() == kHipoAugmentedString) return rowsAS.size(); + return -1; +} +std::string KktMatrix::nla() const { + if (!ptrNE.empty()) return kHipoNormalEqString; + if (!ptrAS.empty()) return kHipoAugmentedString; + return "empty"; +} + } // namespace hipo \ No newline at end of file diff --git a/highs/ipm/hipo/ipm/KktMatrix.h b/highs/ipm/hipo/ipm/KktMatrix.h index c83ea754023..5703b1fb1cf 100644 --- a/highs/ipm/hipo/ipm/KktMatrix.h +++ b/highs/ipm/hipo/ipm/KktMatrix.h @@ -39,6 +39,10 @@ struct KktMatrix { Int buildNEstructure(); Int buildNEvalues(const std::vector& scaling); void freeNEmemory(); + + Int n() const; + Int nz() const; + std::string nla() const; }; } // namespace hipo diff --git a/highs/ipm/hipo/ipm/Solver.cpp b/highs/ipm/hipo/ipm/Solver.cpp index 9a2cc9a66bd..e89ee3686ca 100644 --- a/highs/ipm/hipo/ipm/Solver.cpp +++ b/highs/ipm/hipo/ipm/Solver.cpp @@ -143,7 +143,7 @@ bool Solver::initialise() { } // initialise linear solver - LS_.reset(new FactorHiGHSSolver(*kkt_, options_, model_, regul_, info_, + LS_.reset(new FactorHighsSolver(*kkt_, options_, model_, regul_, info_, it_->data, logger_)); if (!LS_) { info_.status = kStatusError; diff --git a/highs/ipm/hipo/ipm/Solver.h b/highs/ipm/hipo/ipm/Solver.h index 0848667f0df..94e322c2875 100644 --- a/highs/ipm/hipo/ipm/Solver.h +++ b/highs/ipm/hipo/ipm/Solver.h @@ -4,7 +4,7 @@ #include #include "Control.h" -#include "FactorHiGHSSolver.h" +#include "FactorHighsSolver.h" #include "Info.h" #include "Iterate.h" #include "LinearSolver.h" @@ -16,7 +16,7 @@ #include "ipm/hipo/auxiliary/Auxiliary.h" #include "ipm/hipo/auxiliary/IntConfig.h" #include "ipm/hipo/auxiliary/VectorOperations.h" -#include "ipm/hipo/factorhighs/FactorHiGHS.h" +#include "ipm/hipo/factorhighs/FactorHighs.h" #include "ipm/ipx/lp_solver.h" #include "lp_data/HighsCallback.h" #include "lp_data/HighsInfo.h" diff --git a/highs/meson.build b/highs/meson.build index 123afaf937d..dd86d030519 100644 --- a/highs/meson.build +++ b/highs/meson.build @@ -194,7 +194,7 @@ _ipx_srcs = [ _hipo_srcs = [ 'ipm/hipo/ipm/IpmData.cpp', - 'ipm/hipo/ipm/FactorHiGHSSolver.cpp', + 'ipm/hipo/ipm/FactorHighsSolver.cpp', 'ipm/hipo/ipm/Control.cpp', 'ipm/hipo/ipm/Iterate.cpp', 'ipm/hipo/ipm/KktMatrix.cpp', @@ -209,7 +209,8 @@ _hipo_srcs = [ 'ipm/hipo/factorhighs/DenseFactHybrid.cpp', 'ipm/hipo/factorhighs/DenseFactKernel.cpp', 'ipm/hipo/factorhighs/DgemmParallel.cpp', - 'ipm/hipo/factorhighs/FactorHiGHS.cpp', + 'ipm/hipo/factorhighs/FactorHighs_c_api.cpp', + 'ipm/hipo/factorhighs/FactorHighs.cpp', 'ipm/hipo/factorhighs/Factorise.cpp', 'ipm/hipo/factorhighs/FormatHandler.cpp', 'ipm/hipo/factorhighs/HybridHybridFormatHandler.cpp',