Merge branch 'latest' into fix-2460

Author: jajhall

BUILD.bazel

@@ -1,5 +1,7 @@
+load("@bazel_skylib//rules:common_settings.bzl", "bool_flag")
 load("@bazel_skylib//rules:copy_file.bzl", "copy_file")
 load("@rules_cc//cc:defs.bzl", "cc_binary", "cc_library", "cc_test")
+load("@rules_cuda//cuda:defs.bzl", "cuda_library", "requires_cuda")
 
 copy_file(
     name = "highs-config",
@@ -8,12 +10,43 @@ copy_file(
     visibility = ["//visibility:public"],
 )
 
+bool_flag(
+    name = "cupdlp_gpu",
+    build_setting_default = False,
+    visibility = ["//visibility:public"],
+)
+
+config_setting(
+    name = "cupdlp_gpu_enabled",
+    flag_values = {
+        "@rules_cuda//cuda:enable": "True",
+        "@local_cuda//:valid_toolchain_found": "True",
+        ":cupdlp_gpu": "True",
+    },
+)
+
 cc_library(
     name = "config",
     srcs = ["HConfig.h"],
     visibility = ["//visibility:public"],
 )
 
+cuda_library(
+    name = "cupdlp",
+    srcs = glob([
+        "highs/pdlp/cupdlp/cuda/*.cu",
+    ]),
+    hdrs = glob([
+        "highs/pdlp/cupdlp/cuda/*.cuh",
+    ]),
+    defines = ["CUPDLP_GPU"],
+    target_compatible_with = requires_cuda(),
+    deps = [
+        "@local_cuda//:cublas",
+        "@local_cuda//:cusparse",
+    ],
+)
+
 cc_library(
     name = "highs",
     srcs = ["highs/interfaces/highs_c_api.cpp"] + glob([
@@ -49,6 +82,10 @@ cc_library(
             "-Wno-unused-but-set-variable",
         ],
     }),
+    defines = select({
+        ":cupdlp_gpu_enabled": ["CUPDLP_GPU"],
+        "//conditions:default": ["CUPDLP_CPU"],
+    }),
     includes = [
         "extern",
         # "extern/filereaderlp",
@@ -66,7 +103,6 @@ cc_library(
         # "highs/simplex",
         # "highs/test_kkt",
         # "highs/util",
-        "bazel-bin",
     ],
     linkopts = select({
         "@rules_cc//cc/compiler:msvc-cl": ["-DEFAULTLIB:shell32.lib"],
@@ -76,7 +112,10 @@ cc_library(
     deps = [
         "//:config",
         "@zlib",
-    ],
+    ] + select({
+        ":cupdlp_gpu_enabled": [":cupdlp"],
+        "//conditions:default": [],
+    }),
 )
 
 cc_library(

FEATURES.md

@@ -12,4 +12,7 @@ The irreducible infeasibility system (IIS) facility now detects infeasibility du
 
 Prompted by [#2463](https://github.com/ERGO-Code/HiGHS/issues/2463), the HiGHS solution and basis files now match data to any column and row names in the model, only assuming that the data are aligned with column and row indices if there are no names in the model. This requires a new version (v2) of the HiGHS basis file. Basis files from v1 are still read, but deprecated. Now, when writing out a model, basis or solution, column and row names are added to the model - previously they were created temporarily and inconsistentyly on the fly. If the model has existing names, then distinctive names are created to replace any blank names, but names with spaces or duplicate names yield an error status return.
 
+Refactored strong branching to minimize duplicated code
+
 As per [#2487](https://github.com/ERGO-Code/HiGHS/issues/2487), trivial heuristics now run before feasibility jump (FJ), and FJ will use any existing incumbent. FJ will clip any finite variable values in the incumbent to lower and upper bounds, and falls back to the existing logic (lower bound if finite, else upper bound if finite, else 0) for any infinite values in the incumbent.
+

MODULE.bazel

@@ -19,3 +19,14 @@ bazel_dep(
     name = "zlib",
     version = "1.3.1.bcr.5",
 )
+
+bazel_dep(
+    name = "rules_cuda",
+    version = "0.2.5",
+)
+
+local_cuda = use_extension("@rules_cuda//cuda:extensions.bzl", "toolchain")
+local_cuda.local_toolchain(
+    toolkit_path = "",
+)
+use_repo(local_cuda, "local_cuda")

check/TestCAPI.c

@@ -1820,7 +1820,6 @@ void testGetModel() {
   assert(ck_sense == sense);
 
   double* ck_col_cost = (double*)malloc(sizeof(double) * ck_num_col);
-  ;
   double* ck_col_lower = (double*)malloc(sizeof(double) * ck_num_col);
   double* ck_col_upper = (double*)malloc(sizeof(double) * ck_num_col);
   double* ck_row_lower = (double*)malloc(sizeof(double) * ck_num_row);
@@ -2301,6 +2300,112 @@ void testUserObjectiveBoundScaling() {
 
   double scaled_objective_value = Highs_getObjectiveValue(highs);
   assertDoubleValuesEqual("objective_value", unscaled_objective_value, scaled_objective_value);
+void testFixedLp() {
+  // The use of Highs_getFixedLp is illustrated for the MIP
+  //
+  // Min    f  = -3x_0 - 2x_1 - x_2
+  // s.t.          x_0 +  x_1 + x_2 <=  7
+  //              4x_0 + 2x_1 + x_2  = 12
+  //              x_0 >=0; x_1 >= 0; x_2 binary
+
+  const HighsInt num_col = 3;
+  const HighsInt num_row = 2;
+  const HighsInt num_nz = 6;
+  HighsInt a_format = kHighsMatrixFormatColwise;
+  HighsInt sense = kHighsObjSenseMinimize;
+  double offset = 0;
+
+  // Define the column costs, lower bounds and upper bounds
+  double col_cost[3] = {-3.0, -2.0, -1.0};
+  double col_lower[3] = {0.0, 0.0, 0.0};
+  double col_upper[3] = {1.0e30, 1.0e30, 1.0};
+  // Define the row lower bounds and upper bounds
+  double row_lower[2] = {-1.0e30, 12.0};
+  double row_upper[2] = {7.0, 12.0};
+  // Define the constraint matrix column-wise
+  HighsInt a_start[3] = {0, 2, 4};
+  HighsInt a_index[6] = {0, 1, 0, 1, 0, 1};
+  double a_value[6] = {1.0, 4.0, 1.0, 2.0, 1.0, 1.0};
+  HighsInt integrality[3] = {kHighsVarTypeContinuous, kHighsVarTypeContinuous,
+                             kHighsVarTypeInteger};
+
+  void* highs = Highs_create();
+  Highs_setBoolOptionValue(highs, "output_flag", dev_run);
+  Highs_setStringOptionValue(highs, "presolve", "off");
+  HighsInt return_status =
+    Highs_passMip(highs, num_col, num_row, num_nz, a_format, sense, offset,
+		  col_cost, col_lower, col_upper, row_lower, row_upper,
+		  a_start, a_index, a_value, integrality);
+  assert(return_status == kHighsStatusOk);
+  return_status = Highs_run(highs);
+  double mip_objective_function_value;
+  return_status = Highs_getDoubleInfoValue(highs, "objective_function_value",
+                                           &mip_objective_function_value);
+  assert(return_status == kHighsStatusOk);
+
+  double* col_value = (double*)malloc(sizeof(double) * num_col);
+  return_status = Highs_getSolution(highs, col_value, NULL, NULL, NULL);
+  assert(return_status == kHighsStatusOk);
+
+  HighsInt fixed_lp_num_col;
+  HighsInt fixed_lp_num_row;
+  HighsInt fixed_lp_num_nz;
+  HighsInt fixed_lp_sense;
+  double fixed_lp_offset;
+  Highs_getFixedLp(highs, kHighsMatrixFormatColwise, &fixed_lp_num_col, &fixed_lp_num_row,
+		   &fixed_lp_num_nz, &fixed_lp_sense, &fixed_lp_offset, NULL, NULL, NULL, NULL, NULL,
+		   NULL, NULL, NULL);
+
+  assert(fixed_lp_num_col == num_col);
+  assert(fixed_lp_num_row == num_row);
+  assert(fixed_lp_num_nz == num_nz);
+  assert(fixed_lp_sense == sense);
+
+  double* fixed_lp_col_cost = (double*)malloc(sizeof(double) * fixed_lp_num_col);
+  double* fixed_lp_col_lower = (double*)malloc(sizeof(double) * fixed_lp_num_col);
+  double* fixed_lp_col_upper = (double*)malloc(sizeof(double) * fixed_lp_num_col);
+  double* fixed_lp_row_lower = (double*)malloc(sizeof(double) * fixed_lp_num_row);
+  double* fixed_lp_row_upper = (double*)malloc(sizeof(double) * fixed_lp_num_row);
+  HighsInt* fixed_lp_a_start = (HighsInt*)malloc(sizeof(HighsInt) * fixed_lp_num_col);
+  HighsInt* fixed_lp_a_index = (HighsInt*)malloc(sizeof(HighsInt) * fixed_lp_num_nz);
+  double* fixed_lp_a_value = (double*)malloc(sizeof(double) * num_nz);
+
+  // Get the arrays
+  Highs_getFixedLp(highs, kHighsMatrixFormatColwise, &fixed_lp_num_col, &fixed_lp_num_row,
+		   &fixed_lp_num_nz, &fixed_lp_sense, &fixed_lp_offset, fixed_lp_col_cost, fixed_lp_col_lower,
+		   fixed_lp_col_upper, fixed_lp_row_lower, fixed_lp_row_upper, fixed_lp_a_start, fixed_lp_a_index,
+		   fixed_lp_a_value);
+  
+  return_status = Highs_passLp(highs,
+			       fixed_lp_num_col, fixed_lp_num_row, fixed_lp_num_nz,
+			       kHighsMatrixFormatColwise,
+			       fixed_lp_sense, fixed_lp_offset,
+			       fixed_lp_col_cost, fixed_lp_col_lower, fixed_lp_col_upper,
+			       fixed_lp_row_lower, fixed_lp_row_upper,
+			       fixed_lp_a_start, fixed_lp_a_index, fixed_lp_a_value);
+  assert(return_status == kHighsStatusOk);
+
+  return_status = Highs_setSolution(highs, col_value, NULL, NULL, NULL);
+  assert(return_status == kHighsStatusOk);
+
+  return_status = Highs_run(highs);
+  double objective_function_value;
+  return_status = Highs_getDoubleInfoValue(highs, "objective_function_value",
+                                           &objective_function_value);
+  assert(return_status == kHighsStatusOk);
+  assert(objective_function_value == mip_objective_function_value);
+  
+ 
+  free(col_value);
+  free(fixed_lp_col_cost);
+  free(fixed_lp_col_lower);
+  free(fixed_lp_col_upper);
+  free(fixed_lp_row_lower);
+  free(fixed_lp_row_upper);
+  free(fixed_lp_a_start);
+  free(fixed_lp_a_index);
+  free(fixed_lp_a_value);
+  
   Highs_destroy(highs);
 }
 
@@ -2329,6 +2434,7 @@ int main() {
     testDeleteRowResolveWithBasis();
     testIis();
     testUserObjectiveBoundScaling();
+    testFixedLp();
   return 0;
 }
 //  testSetSolution();

check/TestMipSolver.cpp

@@ -1035,3 +1035,157 @@ TEST_CASE("issue-2432", "[highs_test_mip_solver]") {
         "found\n");
   solve(highs, kHighsOffString, require_model_status, optimal_objective);
 }
+
+TEST_CASE("get-fixed-lp", "[highs_test_mip_solver]") {
+  std::string model = "avgas";
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  HighsLp fixed_lp;
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+  REQUIRE(h.getFixedLp(fixed_lp) == HighsStatus::kError);
+
+  model = "flugpl";
+  model_file = std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+  REQUIRE(h.getFixedLp(fixed_lp) == HighsStatus::kError);
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  double mip_optimal_objective = h.getInfo().objective_function_value;
+  HighsSolution solution = h.getSolution();
+
+  // Transform the incumbent MIP into the fixed LP
+  HighsLp mip = h.getLp();
+  std::vector<HighsInt> col_set;
+  std::vector<double> fixed_value;
+  for (HighsInt iCol = 0; iCol < mip.num_col_; iCol++) {
+    if (mip.integrality_[iCol] == HighsVarType::kInteger) {
+      col_set.push_back(iCol);
+      fixed_value.push_back(solution.col_value[iCol]);
+    }
+  }
+  h.clearIntegrality();
+  HighsInt num_set_entries = col_set.size();
+  h.changeColsBounds(num_set_entries, col_set.data(), fixed_value.data(),
+                     fixed_value.data());
+  h.setOptionValue("presolve", kHighsOffString);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+  // In calling changeColsBounds, the incumbent solution was always
+  // cleared, so there was no information from which to construct an
+  // advanced basis. Hence simplex starts from a logical basis and
+  // requires a positive number of iterations (#2556)
+  //
+  // Before code to retain solution if changing the bounds and
+  // solution remains feasible
+  //
+  //  REQUIRE(h.getInfo().simplex_iteration_count > 0);
+  REQUIRE(h.getInfo().simplex_iteration_count == 0);
+
+  // Now, passing the MIP solution, there is information from which to
+  // construct an advanced basis. In the case of flugpl, this is
+  // optimal, so no simplex iterations are required
+  h.clearSolver();
+  h.setSolution(solution);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+  REQUIRE(h.getInfo().simplex_iteration_count == 0);
+
+  // Now re-load the MIP, re-solve, and get the fixed LP
+  REQUIRE(h.passModel(mip) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+
+  REQUIRE(h.getFixedLp(fixed_lp) == HighsStatus::kOk);
+
+  REQUIRE(h.passModel(fixed_lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+
+  // Now run from saved solution (without presolve)
+  h.clearSolver();
+  h.setSolution(solution);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+  REQUIRE(h.getInfo().simplex_iteration_count == 0);
+
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+  // Perturb one of the integer variables for code coverage of
+  // warning: makes fixed LP of flugpl infeasible
+  std::vector<HighsVarType> integrality = h.getLp().integrality_;
+  for (HighsInt iCol = 0; iCol < fixed_lp.num_col_; iCol++) {
+    if (integrality[iCol] != HighsVarType::kContinuous) {
+      solution.col_value[iCol] -= 0.01;
+      break;
+    }
+  }
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  h.setSolution(solution);
+
+  REQUIRE(h.getFixedLp(fixed_lp) == HighsStatus::kWarning);
+
+  REQUIRE(h.passModel(fixed_lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getModelStatus() == HighsModelStatus::kInfeasible);
+
+  h.resetGlobalScheduler(true);
+}
+
+TEST_CASE("get-fixed-lp-semi", "[highs_test_mip_solver]") {
+  HighsLp lp;
+  lp.num_col_ = 4;
+  lp.num_row_ = 2;
+  lp.col_cost_ = {1, 3, 1, 2};
+  lp.col_lower_ = {0, 0, 1, 1};
+  lp.col_upper_ = {1, 1, 3, 5};
+  lp.integrality_ = {HighsVarType::kContinuous, HighsVarType::kInteger,
+                     HighsVarType::kSemiContinuous, HighsVarType::kSemiInteger};
+  lp.row_lower_ = {4, 10};
+  lp.row_upper_ = {kHighsInf, kHighsInf};
+  lp.a_matrix_.start_ = {0, 2, 4, 6, 8};
+  lp.a_matrix_.index_ = {0, 1, 0, 1, 0, 1, 0, 1, 0, 1};
+  lp.a_matrix_.value_ = {1, 1, 1, 2, 1, 3, 1, 4, 5, 1};
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  h.setOptionValue("presolve", kHighsOffString);
+  h.passModel(lp);
+  h.run();
+  double mip_optimal_objective = h.getInfo().objective_function_value;
+  HighsSolution solution = h.getSolution();
+  HighsLp fixed_lp;
+  REQUIRE(h.getFixedLp(fixed_lp) == HighsStatus::kOk);
+
+  REQUIRE(h.passModel(fixed_lp) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+
+  REQUIRE(h.getInfo().objective_function_value == mip_optimal_objective);
+}
+
+TEST_CASE("row-fixed-lp", "[highs_test_mip_solver]") {
+  std::string model = "flugpl";
+  std::string model_file =
+      std::string(HIGHS_DIR) + "/check/instances/" + model + ".mps";
+  Highs h;
+  //  h.setOptionValue("output_flag", dev_run);
+  REQUIRE(h.readModel(model_file) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  double mip_optimal_objective = h.getInfo().objective_function_value;
+  HighsSolution solution = h.getSolution();
+
+  HighsLp lp = h.getLp();
+  h.clearIntegrality();
+  h.changeRowsBounds(0, lp.num_row_ - 1, solution.row_value.data(),
+                     solution.row_value.data());
+  h.setOptionValue("presolve", kHighsOffString);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(h.getInfo().objective_function_value <= mip_optimal_objective);
+
+  h.resetGlobalScheduler(true);
+}