Merged hpm into this branch and resolved minor conflicts

Author: jajhall

CMakeLists.txt

@@ -808,9 +808,6 @@ else(FAST_BUILD)
   if(BUILD_CXX)
     add_subdirectory(app)
     add_subdirectory(examples)
-    if (HIPO)
-      add_subdirectory(highs/ipm/hipo)
-    endif()
   endif()
 
   if(BUILD_TESTING)

check/CMakeLists.txt

@@ -48,6 +48,7 @@ if ((NOT FAST_BUILD OR ALL_TESTS) AND NOT (BUILD_EXTRA_UNIT_ONLY))
       TestHighsRbTree.cpp
       TestHighsSparseMatrix.cpp
       TestHighsVersion.cpp
+      TestHipo.cpp
       TestHSet.cpp
       TestICrash.cpp
       TestIis.cpp

check/TestHipo.cpp

@@ -0,0 +1,74 @@
+#include "HCheckConfig.h"
+#include "Highs.h"
+#include "catch.hpp"
+#include "io/Filereader.h"
+#include "ipm/hipo/ipm/Solver.h"
+#include "ipm/hipo/ipm/Status.h"
+#include "lp_data/HighsCallback.h"
+#include "parallel/HighsParallel.h"
+
+// Example for using HiPO from its C++ interface. The program solves the Netlib
+// problem afiro.
+
+#include <unistd.h>
+
+#include <cmath>
+#include <iostream>
+#include <vector>
+
+const bool dev_run = false;
+
+TEST_CASE("test-hipo-afiro", "[highs_hipo]") {
+  // Test that hipo runs and finds correct solution for afiro
+
+  std::string model = "afiro.mps";
+  const double expected_obj = -464.753;
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue("solver", kHipoString);
+  highs.setOptionValue("timeless_log", kHighsOnString);
+
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/" + model;
+  highs.readModel(filename);
+
+  HighsStatus status = highs.run();
+  REQUIRE(status == HighsStatus::kOk);
+
+  const double actual_obj = highs.getObjectiveValue();
+  REQUIRE(std::abs(actual_obj - expected_obj) < 0.001);
+
+  highs.resetGlobalScheduler(true);
+}
+
+TEST_CASE("test-hipo-deterministic", "[highs_hipo]") {
+  // Test that hipo finds the exact same solution if run twice
+
+  std::string model = "80bau3b.mps";
+
+  Highs highs;
+  highs.setOptionValue("output_flag", dev_run);
+  highs.setOptionValue(kSolverString, kHipoString);
+  highs.setOptionValue(kParallelString, kHighsOnString);
+  highs.setOptionValue(kRunCrossoverString, kHighsOffString);
+
+  std::string filename = std::string(HIGHS_DIR) + "/check/instances/" + model;
+  highs.readModel(filename);
+
+  HighsStatus status = highs.run();
+  REQUIRE(status == HighsStatus::kOk);
+
+  const HighsSolution solution_1 = highs.getSolution();
+
+  highs.run();
+  const HighsSolution solution_2 = highs.getSolution();
+
+  REQUIRE(solution_1.value_valid == solution_2.value_valid);
+  REQUIRE(solution_1.dual_valid == solution_2.dual_valid);
+  REQUIRE(solution_1.col_value == solution_2.col_value);
+  REQUIRE(solution_1.row_value == solution_2.row_value);
+  REQUIRE(solution_1.col_dual == solution_2.col_dual);
+  REQUIRE(solution_1.row_dual == solution_2.row_dual);
+
+  highs.resetGlobalScheduler(true);
+}

cmake/FindHipoDeps.cmake

@@ -3,22 +3,24 @@ option(BLAS_ROOT "Root directory of BLAS or OpenBLAS" "")
 message(STATUS "BLAS_ROOT is " ${BLAS_ROOT})
 
 if(NOT APPLE)
-    find_library(BLAS_LIB 
-        NAMES blas HINTS 
-        "${BLAS_ROOT}/lib")
+    find_library(OPENBLAS_LIB 
+        NAMES openblas 
+        HINTS "${BLAS_ROOT}/lib")
 
-    if(BLAS_LIB)
-        message("Found BLAS library at ${BLAS_LIB}")
-    else(BLAS_LIB)
-        find_library(OPENBLAS_LIB 
-            NAMES openblas 
-            HINTS "${BLAS_ROOT}/lib")
-        if(NOT OPENBLAS_LIB)
-            message(FATAL_ERROR "No BLAS library found")
-        endif(NOT OPENBLAS_LIB)
+    if(OPENBLAS_LIB)
         message("Found OpenBLAS library at ${OPENBLAS_LIB}")
-    endif(BLAS_LIB)
-    
+
+    else(OPENBLAS_LIB)
+        find_library(BLAS_LIB 
+            NAMES blas HINTS 
+            "${BLAS_ROOT}/lib")
+
+        if(NOT BLAS_LIB)
+            message(FATAL_ERROR "No BLAS library found")
+        endif(NOT BLAS_LIB)
+        message("Found BLAS library at ${BLAS_LIB}")
+    endif(OPENBLAS_LIB)
+
 endif(NOT APPLE)
 
 # METIS

highs/CMakeLists.txt

@@ -188,12 +188,13 @@ else()
 
     if(APPLE)
         target_link_libraries(highs "-framework Accelerate")
-        target_compile_definitions(highs PRIVATE FRAMEWORK_ACCELERATE)
+        target_compile_definitions(highs PRIVATE HIPO_USES_APPLE_BLAS)
     else(APPLE)
         if(BLAS_LIB)
             target_link_libraries(highs "${BLAS_LIB}")
         elseif(OPENBLAS_LIB)
             target_link_libraries(highs "${OPENBLAS_LIB}")
+            target_compile_definitions(highs PRIVATE HIPO_USES_OPENBLAS)
         else()
             message(FATAL_ERROR "No BLAS library available")
         endif(BLAS_LIB)

highs/ipm/IpxWrapper.cpp

@@ -402,6 +402,13 @@ HighsStatus solveLpHipo(HighsLpSolverObject& solver_object) {
                      solver_object.highs_info_, solver_object.callback_);
 }
 
+#ifdef HIPO_USES_OPENBLAS
+// function to set number of threads of openblas
+extern "C" {
+void openblas_set_num_threads(int num_threads);
+}
+#endif
+
 HighsStatus solveLpHipo(const HighsOptions& options, HighsTimer& timer,
                         const HighsLp& lp, HighsBasis& highs_basis,
                         HighsSolution& highs_solution,
@@ -439,6 +446,11 @@ HighsStatus solveLpHipo(const HighsOptions& options, HighsTimer& timer,
   // Indicate that no imprecise solution has (yet) been found
   resetModelStatusAndHighsInfo(model_status, highs_info);
 
+#ifdef HIPO_USES_OPENBLAS
+  // force openblas to run in serial, for determinism and better performance
+  openblas_set_num_threads(1);
+#endif
+
   // Create solver instance
   hipo::Solver hipo{};
 
@@ -499,13 +511,22 @@ HighsStatus solveLpHipo(const HighsOptions& options, HighsTimer& timer,
     hipo_options.parallel = hipo::kOptionParallelChoose;
   }
 
-  // Highs does not have an option to select NE/AS approach.
-  // For now use choose, but an option should be added for the user to choose.
-  hipo_options.nla = hipo::kOptionNlaChoose;
+  // Parse hipo_system option
+  if (options.hipo_system == kHipoAugmentedString) {
+    hipo_options.nla = hipo::kOptionNlaAugmented;
+  } else if (options.hipo_system == kHipoNormalEqString) {
+    hipo_options.nla = hipo::kOptionNlaNormEq;
+  } else if (options.hipo_system == kHighsChooseString) {
+    hipo_options.nla = hipo::kOptionNlaChoose;
+  } else {
+    highsLogUser(options.log_options, HighsLogType::kError,
+                 "Unknown value of option hipo_system\n");
+    model_status = HighsModelStatus::kSolveError;
+    return HighsStatus::kError;
+  }
 
   // ===========================================================================
   // TO DO
-  // - add options for NE/AS
   // - consider adding options for parallel tree/node
   // - block size for dense factorisation can have large impact on performance
   //   and depends on the specific architecture. It may be worth exposing it to

highs/ipm/hipo/CMakeLists.txt

@@ -588,7 +588,7 @@ void FactorHiGHSSolver::setParallel(Options& options) {
       log_stream << "Full requested\n";
       break;
     case kOptionParallelChoose: {
-#ifdef FRAMEWORK_ACCELERATE
+#ifdef HIPO_USES_APPLE_BLAS
       // Blas on Apple do not work well with parallel_node, but parallel_tree
       // seems to always be beneficial.
       parallel_node = false;

highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp

@@ -293,6 +293,11 @@ const string kReadSolutionFileString = "read_solution_file";
 // String for HiGHS log file option
 const string kLogFileString = "log_file";
 
+// Strings for HiPO system option
+const string kHipoSystemString = "hipo_system";
+const string kHipoAugmentedString = "augmented";
+const string kHipoNormalEqString = "normaleq";
+
 // Strings for MIP LP/IPM options
 const string kMipLpSolverString = "mip_lp_solver";
 const string kMipIpmSolverString = "mip_ipm_solver";
@@ -359,6 +364,7 @@ struct HighsOptionsStruct {
   // Options for IPM solver
   double ipm_optimality_tolerance;
   HighsInt ipm_iteration_limit;
+  std::string hipo_system;
 
   // Options for PDLP solver
   bool pdlp_scaling;
@@ -527,6 +533,7 @@ struct HighsOptionsStruct {
         timeless_log(false),
         ipm_optimality_tolerance(0.0),
         ipm_iteration_limit(0),
+        hipo_system(""),
         pdlp_scaling(false),
         pdlp_iteration_limit(0),
         pdlp_e_restart_method(0),
@@ -1243,6 +1250,12 @@ class HighsOptions : public HighsOptionsStruct {
         &ipm_iteration_limit, 0, kHighsIInf, kHighsIInf);
     records.push_back(record_int);
 
+    record_string = new OptionRecordString(
+        kHipoSystemString,
+        "HiPO Newton system option: \"augmented\", \"normaleq\" or \"choose\".",
+        advanced, &hipo_system, kHighsChooseString);
+    records.push_back(record_string);
+
     record_bool = new OptionRecordBool(
         "pdlp_scaling", "Scaling option for PDLP solver: Default = true",
         advanced, &pdlp_scaling, true);