Added testing of passing square Hessian via passModel and passHessian

check/TestQpSolver.cpp

@@ -1318,6 +1318,50 @@ TEST_CASE("issue-2821", "[qpsolver]") {
   h.resetGlobalScheduler(true);
 }
 
+TEST_CASE("pass-square-hessian", "[qpsolver]") {
+  Highs h;
+  h.setOptionValue("output_flag", dev_run);
+  const HighsInfo& info = h.getInfo();
+  // Set up the same constraints as 2821, but use a Hessian with more
+  // entries and varied values
+  HighsModel model;
+  HighsLp& lp = model.lp_;
+  HighsHessian& hessian = model.hessian_;
+  lp.num_col_ = 5;
+  lp.num_row_ = 3;
+  lp.col_cost_ = {0, -4, -4, -2, -2};
+  lp.col_lower_ = {0, 0, 0, 0, 0};
+  lp.col_upper_ = {inf, inf, inf, inf, inf};
+  lp.row_lower_ = {4, 0, 0};
+  lp.row_upper_ = {4, 0, 0};
+  lp.a_matrix_.start_ = {0, 1, 3, 4, 5, 7};
+  lp.a_matrix_.index_ = {0, 0, 2, 1, 1, 1, 2};
+  lp.a_matrix_.value_ = {1, 3, 1, 1, 1, -2, -1};
+  hessian.dim_ = lp.num_col_;
+  hessian.start_ = {0, 3, 6, 9, 11, 12};
+  hessian.index_ = {0, 1, 3, 1, 2, 3, 2, 3, 4, 3, 4, 4};
+  hessian.value_ = {4, -1, 1, 8, -2, 3, 9, -3, 2, 6, -4, 5};
+  const double optimal_objective_value = -1.9875776398e-01;
+  REQUIRE(h.passModel(model) == HighsStatus::kOk);
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(okValueDifference(info.objective_function_value,
+                            optimal_objective_value));
+  HighsHessian square_hessian = h.getModel().hessian_.toSquare();
+  REQUIRE(h.passHessian(square_hessian) == HighsStatus::kOk);
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(okValueDifference(info.objective_function_value,
+                            optimal_objective_value));
+  model.hessian_ = square_hessian;
+  REQUIRE(h.passModel(model) == HighsStatus::kOk);
+
+  REQUIRE(h.run() == HighsStatus::kOk);
+  REQUIRE(okValueDifference(info.objective_function_value,
+                            optimal_objective_value));
+
+  h.resetGlobalScheduler(true);
+}
+
 /*
 // This case causes a segfault with the meson build, but not with
 // cmake. Might be indicative of the bug causing the degeneracy

highs/model/HighsHessian.cpp

@@ -143,12 +143,12 @@ HighsInt HighsHessian::numNz() const {
   return this->start_[this->dim_];
 }
 
-void HighsHessian::print() const {
+void HighsHessian::print(const std::string& message) const {
   HighsInt num_nz = this->numNz();
   printf("%s Hessian of dimension %" HIGHSINT_FORMAT " and %" HIGHSINT_FORMAT
-         " entries\n",
+         " entries: %s\n",
          this->format_ == HessianFormat::kTriangular ? "Triangular" : "Square",
-         dim_, num_nz);
+         dim_, num_nz, message.c_str());
   printf("Start; Index; Value of sizes %d; %d; %d\n", (int)this->start_.size(),
          (int)this->index_.size(), (int)this->value_.size());
   if (dim_ <= 0) return;
@@ -256,8 +256,64 @@ bool HighsHessian::isDiagonal() const {
   }
   return true;
 }
+
 double HighsHessian::diag(HighsInt i) const {
   assert(i < dim_);
   assert(index_[start_[i]] == i);
   return value_[start_[i]];
 }
+
+HighsHessian HighsHessian::toSquare() const {
+  if (this->format_ == HessianFormat::kSquare) return *this;
+  assert(this->format_ == HessianFormat::kTriangular);
+  std::vector<HighsInt> iwork(this->dim_, 0);
+  for (HighsInt iCol = 0; iCol < this->dim_; iCol++) {
+    HighsInt iEl = this->start_[iCol];
+    HighsInt iRow = this->index_[iEl];
+    assert(iRow == iCol);
+    iEl++;
+    for (; iEl < this->start_[iCol + 1]; iEl++) {
+      HighsInt iRow = this->index_[iEl];
+      assert(iRow > iCol);
+      iwork[iRow]++;
+    }
+  }
+  HighsInt triangular_hessian_off_diagonal = this->numNz() - this->dim_;
+  HighsHessian square_hessian;
+  square_hessian.format_ = HessianFormat::kSquare;
+  square_hessian.dim_ = this->dim_;
+  square_hessian.start_[0] = 0;
+  for (HighsInt iCol = 0; iCol < this->dim_; iCol++) {
+    HighsInt square_hessian_col_nnz =
+        iwork[iCol] + this->start_[iCol + 1] - this->start_[iCol];
+    square_hessian.start_.push_back(square_hessian.start_[iCol] +
+                                    square_hessian_col_nnz);
+    iwork[iCol] = square_hessian.start_[iCol] + 1;
+  }
+  HighsInt square_hessian_nnz = square_hessian.numNz();
+  assert(square_hessian_nnz ==
+         this->dim_ + 2 * triangular_hessian_off_diagonal);
+  square_hessian.index_.resize(square_hessian_nnz);
+  square_hessian.value_.resize(square_hessian_nnz);
+  for (HighsInt iCol = 0; iCol < this->dim_; iCol++) {
+    HighsInt iEl = this->start_[iCol];
+    HighsInt iRow = this->index_[iEl];
+    HighsInt square_iEl = square_hessian.start_[iCol];
+    square_hessian.index_[square_iEl] = iCol;
+    square_hessian.value_[square_iEl] = this->value_[iEl];
+    iEl++;
+    for (; iEl < this->start_[iCol + 1]; iEl++) {
+      HighsInt iRow = this->index_[iEl];
+      double value = this->value_[iEl];
+      HighsInt square_iEl = iwork[iCol];
+      square_hessian.index_[square_iEl] = iRow;
+      square_hessian.value_[square_iEl] = value;
+      iwork[iCol]++;
+      square_iEl = iwork[iRow];
+      square_hessian.index_[square_iEl] = iCol;
+      square_hessian.value_[square_iEl] = value;
+      iwork[iRow]++;
+    }
+  }
+  return square_hessian;
+}

highs/model/HighsHessian.h

@@ -47,8 +47,8 @@ class HighsHessian {
   bool empty() const;
   bool isDiagonal() const;
   double diag(HighsInt i) const;
-
-  void print() const;
+  HighsHessian toSquare() const;
+  void print(const std::string& message = "") const;
 };
 
 #endif