Add C interface to HiPO linear solver

.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

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();
 }

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

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

highs/ipm/hipo/auxiliary/Auxiliary.cpp

@@ -26,21 +26,19 @@ void subtreeSize(const std::vector<Int>& parent, std::vector<Int>& sizes) {
   }
 }
 
-void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-               std::vector<Int>& ptrT, std::vector<Int>& 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<Int> 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<Int>& ptr, const std::vector<Int>& rows,
 }
 
 void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
-               const std::vector<double>& val, std::vector<Int>& ptrT,
-               std::vector<Int>& rowsT, std::vector<double>& valT) {
-  // Compute the transpose of the matrix and return it in rowsT, ptrT and valT
+               std::vector<Int>& ptrT, std::vector<Int>& 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<Int> 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<Int>& ptr, const std::vector<Int>& rows,
   }
 }
 
+void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
+               const std::vector<double>& val, std::vector<Int>& ptrT,
+               std::vector<Int>& rowsT, std::vector<double>& valT) {
+  transpose(ptr.size() - 1, rows.size(), ptr.data(), rows.data(), val.data(),
+            ptrT.data(), rowsT.data(), valT.data());
+}
+
 void permuteSym(const std::vector<Int>& iperm, std::vector<Int>& ptr,
                 std::vector<Int>& rows, std::vector<double>& 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<Int> work(n, 0);
@@ -117,7 +125,7 @@ void permuteSym(const std::vector<Int>& iperm, std::vector<Int>& 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<Int> new_rows(new_ptr.back());
   std::vector<double> new_val;
@@ -300,16 +308,15 @@ Int maxDepthTree(const std::vector<Int>& parent) {
   return max_depth;
 }
 
-void fullFromLower(const std::vector<Int>& ptrL, const std::vector<Int>& rowsL,
+void fullFromLower(Int n, Int nz, const Int* ptrL, const Int* rowsL,
                    std::vector<Int>& ptrF, std::vector<Int>& rowsF) {
   // Given a sparse matrix in lower triangular format, build the same matrix in
   // full format, without diagonal entries.
 
-  std::vector<Int> rowsU(rowsL.size());
-  std::vector<Int> ptrU(ptrL.size());
-  transpose(ptrL, rowsL, ptrU, rowsU);
+  std::vector<Int> rowsU(nz);
+  std::vector<Int> ptrU(n + 1);
+  transpose(n, nz, ptrL, rowsL, ptrU.data(), rowsU.data());
 
-  const Int n = ptrL.size() - 1;
   std::vector<Int> 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<Int>& ptrL, const std::vector<Int>& 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) {

highs/ipm/hipo/auxiliary/Auxiliary.h

@@ -14,8 +14,12 @@ namespace hipo {
 
 void inversePerm(const std::vector<Int>& perm, std::vector<Int>& iperm);
 void subtreeSize(const std::vector<Int>& parent, std::vector<Int>& sizes);
+void transpose(Int n, Int nz, const Int* ptr, const Int* rows, Int* ptrT,
+               Int* rowsT);
 void transpose(const std::vector<Int>& ptr, const std::vector<Int>& rows,
                std::vector<Int>& ptrT, std::vector<Int>& 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<Int>& ptr, const std::vector<Int>& rows,
                const std::vector<double>& val, std::vector<Int>& ptrT,
                std::vector<Int>& rowsT, std::vector<double>& valT);
@@ -32,19 +36,18 @@ void processEdge(Int j, Int i, const std::vector<Int>& first,
 Int64 getDiagStart(Int n, Int k, Int nb, Int n_blocks,
                    std::vector<Int64>& start, bool triang = false);
 Int maxDepthTree(const std::vector<Int>& parent);
-void fullFromLower(const std::vector<Int>& ptrL, const std::vector<Int>& rowsL,
+void fullFromLower(Int n, Int nz, const Int* ptrL, const Int* rowsL,
                    std::vector<Int>& ptrF, std::vector<Int>& rowsF);
 double snFlops(double size, double clique_size);
 double snSpops(double clique_size);
 
 template <typename T>
-void counts2Ptr(std::vector<T>& ptr, std::vector<T>& 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,13 +56,29 @@ void counts2Ptr(std::vector<T>& ptr, std::vector<T>& w) {
   ptr[n] = temp_nz;
 }
 
+template <typename T>
+void permuteVector(T* v, const std::vector<Int>& perm) {
+  // Permute vector v according to permutation perm.
+  const Int n = perm.size();
+  std::vector<T> temp_v(v, v + n);
+  for (Int i = 0; i < n; ++i) v[i] = temp_v[perm[i]];
+}
+
 template <typename T>
 void permuteVector(std::vector<T>& v, const std::vector<Int>& perm) {
   // Permute vector v according to permutation perm.
   std::vector<T> temp_v(v);
   for (Int i = 0; i < static_cast<Int>(v.size()); ++i) v[i] = temp_v[perm[i]];
 }
 
+template <typename T>
+void permuteVectorInverse(T* v, const std::vector<Int>& iperm) {
+  // Permute vector v according to inverse permutation iperm.
+  const Int n = iperm.size();
+  std::vector<T> temp_v(v, v + n);
+  for (Int i = 0; i < n; ++i) v[iperm[i]] = temp_v[i];
+}
+
 template <typename T>
 void permuteVectorInverse(std::vector<T>& v, const std::vector<Int>& iperm) {
   // Permute vector v according to inverse permutation iperm.

highs/ipm/hipo/auxiliary/Logger.cpp

@@ -1,7 +1,11 @@
 #include "Logger.h"
 
+#include <stdarg.h>
+
 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