Remove redundant copy operations in FactorHiGHSSolver

Add Zed editor caches to .gitignore
Fix a compiler warning about potential infinite loops

Author: aasivas

.gitignore

@@ -33,6 +33,10 @@ local.properties
 # PDT-specific
 .buildpath
 
+#################
+## Zed Editor
+#################
+.cache
 
 #################
 ## Visual Studio
@@ -292,4 +296,4 @@ highs-problem-set
 # Generated HiGHS options file
 Options.md
 
-build-latest
+build-latest

highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp

@@ -1,7 +1,5 @@
 #include "FactorHiGHSSolver.h"
 
-#include <limits>
-
 #include "Status.h"
 #include "ipm/hipo/auxiliary/Auxiliary.h"
 #include "ipm/hipo/auxiliary/Log.h"
@@ -92,9 +90,9 @@ Int FactorHiGHSSolver::buildNEvalues(const HighsSparseMatrix& A,
                                      const std::vector<double>& scaling) {
   // given the NE structure already computed, fill in the NE values
 
-  assert(!ptrNE_.empty() && !rowsNE_.empty());
+  assert(!ptrLower_.empty() && !rowsLower_.empty());
 
-  valNE_.assign(rowsNE_.size(), 0.0);
+  valLower_.assign(rowsLower_.size(), 0.0);
 
   std::vector<double> work(A.num_row_, 0.0);
 
@@ -127,9 +125,9 @@ Int FactorHiGHSSolver::buildNEvalues(const HighsSparseMatrix& A,
     // intersection of row with rows below finished.
 
     // read from work, using indices of column "row" of AAt
-    for (Int el = ptrNE_[row]; el < ptrNE_[row + 1]; ++el) {
-      Int index = rowsNE_[el];
-      valNE_[el] = work[index];
+    for (Int el = ptrLower_[row]; el < ptrLower_[row + 1]; ++el) {
+      Int index = rowsLower_[el];
+      valLower_[el] = work[index];
       work[index] = 0.0;
     }
   }
@@ -170,11 +168,11 @@ Int FactorHiGHSSolver::buildNEstructure(const HighsSparseMatrix& A,
     }
   }
 
-  ptrNE_.clear();
-  rowsNE_.clear();
+  ptrLower_.clear();
+  rowsLower_.clear();
 
   // ptr is allocated its exact size
-  ptrNE_.resize(A.num_row_ + 1, 0);
+  ptrLower_.resize(A.num_row_ + 1, 0);
 
   // keep track if given entry is nonzero, in column considered
   std::vector<bool> is_nz(A.num_row_, false);
@@ -211,74 +209,101 @@ Int FactorHiGHSSolver::buildNEstructure(const HighsSparseMatrix& A,
     // intersection of row with rows below finished.
 
     // if the total number of nonzeros exceeds the maximum, return error.
-    if ((Int64)ptrNE_[row] + nz_in_col >= nz_limit) return kStatusOverflow;
+    if ((Int64)ptrLower_[row] + nz_in_col >= nz_limit) return kStatusOverflow;
 
     // update pointers
-    ptrNE_[row + 1] = ptrNE_[row] + nz_in_col;
+    ptrLower_[row + 1] = ptrLower_[row] + nz_in_col;
 
     // now assign indices
     for (Int i = 0; i < nz_in_col; ++i) {
       Int index = temp_index[i];
       // push_back is better then reserve, because the final length is not known
-      rowsNE_.push_back(index);
+      rowsLower_.push_back(index);
       is_nz[index] = false;
     }
   }
 
   return kStatusOk;
 }
 
-Int FactorHiGHSSolver::factorAS(const HighsSparseMatrix& A,
-                                const std::vector<double>& scaling) {
-  // only execute factorisation if it has not been done yet
-  assert(!this->valid_);
-
-  Clock clock;
-
-  // initialise
-  std::vector<Int> ptrLower;
-  std::vector<Int> rowsLower;
-  std::vector<double> valLower;
-
+Int FactorHiGHSSolver::buildASstructure(const HighsSparseMatrix& A,
+                                        const std::vector<double>& scaling) {
   Int nA = A.num_col_;
   Int mA = A.num_row_;
   Int nzA = A.numNz();
 
-  ptrLower.resize(nA + mA + 1);
-  rowsLower.resize(nA + nzA + mA);
-  valLower.resize(nA + nzA + mA);
+  ptrLower_.resize(nA + mA + 1);
+  rowsLower_.resize(nA + nzA + mA);
+  valLower_.resize(nA + nzA + mA);
 
   // build lower triangle
   Int next = 0;
 
   for (Int i = 0; i < nA; ++i) {
     // diagonal element
-    rowsLower[next] = i;
-    valLower[next++] = -scaling[i];
+    rowsLower_[next] = i;
+    valLower_[next] = -scaling[i];
+    next++;
 
     // column of A
     for (Int el = A.start_[i]; el < A.start_[i + 1]; ++el) {
-      rowsLower[next] = A.index_[el] + nA;
-      valLower[next++] = A.value_[el];
+      rowsLower_[next] = A.index_[el] + nA;
+      valLower_[next] = A.value_[el];
+      next++;
     }
 
-    ptrLower[i + 1] = next;
+    ptrLower_[i + 1] = next;
   }
 
   // 2,2 block
   for (Int i = 0; i < mA; ++i) {
-    rowsLower[next] = nA + i;
-    valLower[next++] = 0.0;
-    ptrLower[nA + i + 1] = ptrLower[nA + i] + 1;
+    rowsLower_[next] = nA + i;
+    valLower_[next++] = 0.0;
+    ptrLower_[nA + i + 1] = ptrLower_[nA + i] + 1;
+  }
+  return 0;
+}
+
+Int FactorHiGHSSolver::updateASdiagonal(const HighsSparseMatrix& A,
+                                      const std::vector<double>& scaling) {
+  Int nA = A.num_col_;
+
+  // build lower triangle
+  Int next = 0;
+
+  for (Int i = 0; i < nA; ++i) {
+    // diagonal element
+    valLower_[next] = -scaling[i];
+    next++;
+
+    // move the pointer by the number of nonzeros in the i-th column of A
+    next += A.start_[i + 1] - A.start_[i];
+  }
+  return 0;
+}
+
+Int FactorHiGHSSolver::factorAS(const HighsSparseMatrix& A,
+                                const std::vector<double>& scaling) {
+  // only execute factorisation if it has not been done yet
+  assert(!this->valid_);
+
+  Clock clock;
+
+  if (amat_set_ == true) {
+    std::ignore = updateASdiagonal(A, scaling);
+  } else {
+    std::ignore = buildASstructure(A, scaling);
+    amat_set_ = true;
   }
+
   if (info_) info_->matrix_time += clock.stop();
 
   // set static regularisation, since it may have changed
   FH_.setRegularisation(regul_.primal, regul_.dual);
 
   // factorise matrix
   clock.start();
-  if (FH_.factorise(S_, rowsLower, ptrLower, valLower))
+  if (FH_.factorise(S_, rowsLower_, ptrLower_, valLower_))
     return kStatusErrorFactorise;
   if (info_) {
     info_->factor_time += clock.stop();
@@ -305,11 +330,8 @@ Int FactorHiGHSSolver::factorNE(const HighsSparseMatrix& A,
 
   // factorise
   clock.start();
-  // make copies of structure, because factorise object will take ownership and
-  // modify them.
-  std::vector<Int> ptrNE(ptrNE_);
-  std::vector<Int> rowsNE(rowsNE_);
-  if (FH_.factorise(S_, rowsNE, ptrNE, valNE_)) return kStatusErrorFactorise;
+  if (FH_.factorise(S_, rowsLower_, ptrLower_, valLower_))
+    return kStatusErrorFactorise;
   if (info_) {
     info_->factor_time += clock.stop();
     info_->factor_number++;
@@ -385,9 +407,7 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) {
   log_.printDevInfo("Building AS structure\n");
 
   Clock clock;
-  std::vector<Int> ptrLower;
-  std::vector<Int> rowsLower;
-  if (Int status = getASstructure(model_.A(), ptrLower, rowsLower))
+  if (Int status = getASstructure(model_.A(), ptrLower_, rowsLower_))
     return status;
   if (info_) info_->matrix_structure_time = clock.stop();
 
@@ -399,7 +419,7 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) {
   log_.printDevInfo("Performing AS analyse phase\n");
 
   clock.start();
-  Int status = FH_.analyse(S, rowsLower, ptrLower, pivot_signs);
+  Int status = FH_.analyse(S, rowsLower_, ptrLower_, pivot_signs);
   if (info_) info_->analyse_AS_time = clock.stop();
 
   if (status && log_.debug(2)) {
@@ -410,17 +430,6 @@ Int FactorHiGHSSolver::analyseAS(Symbolic& S) {
   return status ? kStatusErrorAnalyse : kStatusOk;
 }
 
-void FactorHiGHSSolver::freeNEmemory() {
-  // Swap NE data structures with empty vectors, to guarantee that memory is
-  // freed.
-
-  std::vector<Int>().swap(ptrNE_);
-  std::vector<Int>().swap(rowsNE_);
-  std::vector<Int>().swap(ptrA_rw_);
-  std::vector<Int>().swap(idxA_rw_);
-  std::vector<Int>().swap(corr_A_);
-}
-
 Int FactorHiGHSSolver::analyseNE(Symbolic& S, Int64 nz_limit) {
   // Perform analyse phase of augmented system and return symbolic factorisation
   // in object S and the status. If building the matrix failed, the status is
@@ -438,7 +447,7 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S, Int64 nz_limit) {
   log_.printDevInfo("Performing NE analyse phase\n");
 
   clock.start();
-  Int status = FH_.analyse(S, rowsNE_, ptrNE_, pivot_signs);
+  Int status = FH_.analyse(S, rowsLower_, ptrLower_, pivot_signs);
   if (info_) info_->analyse_NE_time = clock.stop();
 
   if (status && log_.debug(2)) {
@@ -548,7 +557,6 @@ Int FactorHiGHSSolver::chooseNla() {
   if (status == kStatusOk) {
     if (options_.nla == kOptionNlaAugmented) {
       S_ = std::move(symb_AS);
-      freeNEmemory();
     } else {
       S_ = std::move(symb_NE);
     }
@@ -681,4 +689,4 @@ void FactorHiGHSSolver::setParallel() {
   S_.setParallel(parallel_tree, parallel_node);
 }
 
-}  // namespace hipo
+}  // namespace hipo

highs/ipm/hipo/ipm/FactorHiGHSSolver.h

@@ -1,8 +1,6 @@
 #ifndef HIPO_FACTORHIGHS_SOLVER_H
 #define HIPO_FACTORHIGHS_SOLVER_H
 
-#include <algorithm>
-
 #include "Info.h"
 #include "IpmData.h"
 #include "LinearSolver.h"
@@ -20,11 +18,12 @@ class FactorHiGHSSolver : public LinearSolver {
   Symbolic S_;
 
   // normal equations data
-  std::vector<Int> ptrNE_;
-  std::vector<Int> rowsNE_;
-  std::vector<double> valNE_;
-  std::vector<Int> ptrA_rw_, idxA_rw_;
-  std::vector<Int> corr_A_;
+  std::vector<Int> ptrLower_;
+  std::vector<Int> rowsLower_;
+  std::vector<double> valLower_;
+  std::vector<Int> ptrA_rw_, idxA_rw_; // Symbolic row-wise (CSR) copy of A
+  std::vector<Int> corr_A_; // ?
+  bool amat_set_ = false;
 
   const Regularisation& regul_;
 
@@ -42,7 +41,11 @@ class FactorHiGHSSolver : public LinearSolver {
   Int buildNEstructure(const HighsSparseMatrix& A, Int64 nz_limit = kHighsIInf);
   Int buildNEvalues(const HighsSparseMatrix& A,
                     const std::vector<double>& scaling);
-  void freeNEmemory();
+
+  Int buildASstructure(const HighsSparseMatrix& A,
+                       const std::vector<double>& scaling);
+  Int updateASdiagonal(const HighsSparseMatrix& A,
+      const std::vector<double>& scaling);
 
   Int analyseAS(Symbolic& S);
   Int analyseNE(Symbolic& S, Int64 nz_limit = kHighsIInf);

highs/ipm/hipo/ipm/Solver.cpp

@@ -2,6 +2,7 @@
 
 #include <cassert>
 #include <cmath>
+#include <cstdint>
 #include <iostream>
 
 #include "ipm/hipo/auxiliary/AutoDetect.h"
@@ -417,7 +418,7 @@ double Solver::stepToBoundary(const std::vector<double>& x,
   double alpha = 1.0;
   Int bl = -1;
 
-  for (Int i = 0; i < x.size(); ++i) {
+  for (uint64_t i = 0; i < x.size(); ++i) {
     if ((lo && model_.hasLb(i)) || (!lo && model_.hasUb(i))) {
       double c = (cor ? (*cor)[i] * weight : 0.0);
       if (x[i] + alpha * (dx[i] + c) < 0.0) {