NE structure now uses linked lists to access rows

filikat · filikat · commit a98088454274 · 2025-09-18T15:45:13.000+01:00
diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp b/highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp
@@ -72,80 +72,92 @@ Int FactorHiGHSSolver::setup() {
   return kStatusOk;
 }
 
-Int FactorHiGHSSolver::buildNEstructureDense(const HighsSparseMatrix& A,
-                                             int64_t nz_limit) {
-  // Build lower triangular structure of AAt.
-  // This approach should be advantageous if AAt is expected to be relatively
-  // dense.
-
-  // create row-wise copy of the matrix
-  AT_ = A;
-  AT_.ensureRowwise();
-
-  ptrNE_.clear();
-  rowsNE_.clear();
+Int FactorHiGHSSolver::buildNEvalues(const HighsSparseMatrix& A,
+                                     const std::vector<double>& scaling) {
+  // given the NE structure already computed, fill in the NE values
 
-  // ptr is allocated its exact size
-  ptrNE_.resize(A.num_row_ + 1, 0);
+  assert(!ptrNE_.empty() && !rowsNE_.empty());
 
-  // temporary dense vector
-  std::vector<Int> work(A.num_row_);
+  valNE_.assign(rowsNE_.size(), 0.0);
 
-  int64_t AAt_nz = 0;
+  std::vector<double> work(A.num_row_, 0.0);
 
   for (Int row = 0; row < A.num_row_; ++row) {
-    // work contains information about column "row".
-    // if there is a 1 in position pos, then AAt has nonzero in entry (pos,row).
-    // only entries below entry "row" (the diagonal) are used.
-
     // go along the entries of the row, and then down each column.
     // this builds the lower triangular part of the row-th column of AAt.
     for (Int rowEl = AT_.start_[row]; rowEl < AT_.start_[row + 1]; ++rowEl) {
       Int col = AT_.index_[rowEl];
 
+      const double theta =
+          scaling.empty() ? 1.0 : 1.0 / (scaling[col] + regul_.primal);
+
+      const double row_value = theta * AT_.value_[rowEl];
+
       // for each nonzero in the row, go down corresponding column
       for (Int colEl = A.start_[col]; colEl < A.start_[col + 1]; ++colEl) {
         Int row2 = A.index_[colEl];
 
         // skip when row2 is above row
         if (row2 < row) continue;
 
-        // save information that there is nonzero in position (row2,row).
-        work[row2] = 1;
+        // compute and accumulate value
+        double value = row_value * A.value_[colEl];
+        work[row2] += value;
       }
     }
     // intersection of row with rows below finished.
-    // now work contains the sparsity pattern of AAt(row:end,row).
 
-    // now assign indices
-    Int col_nz = 0;
-    for (Int i = row; i < work.size(); ++i) {
-      if (work[i]) {
-        if (AAt_nz + 1 >= nz_limit) return kStatusOoM;
-
-        rowsNE_.push_back(i);
-        work[i] = 0;
-        ++AAt_nz;
-        ++col_nz;
-      }
+    // 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];
+      work[index] = 0.0;
     }
-
-    // update pointers
-    ptrNE_[row + 1] = ptrNE_[row] + col_nz;
   }
 
   return kStatusOk;
 }
 
-Int FactorHiGHSSolver::buildNEstructureSparse(const HighsSparseMatrix& A,
-                                              int64_t nz_limit) {
+Int FactorHiGHSSolver::buildNEstructure(const HighsSparseMatrix& A,
+                                        int64_t nz_limit) {
   // Build lower triangular structure of AAt.
-  // This approach should be advantageous if AAt is expected to be very sparse.
+  // This approach uses a column-wise copy of A and a collection of linked lists
+  // to access the rows
 
-  // create row-wise copy of the matrix
+  // create row-wise copy of the matrix (to be removed later!!!!!!!!)
   AT_ = A;
   AT_.ensureRowwise();
 
+  // NB: A must have sorted columns for this to work
+
+  // colNE stores the index of the column of each entry.
+  // It's useless when accessing the matrix by columns, but useful when doing it
+  // by rows.
+  colNE_.resize(A.numNz());
+
+  // head and next are a collection of linked lists to access the matrix by
+  // rows.
+  nextNE_.resize(A.numNz(), -1);
+  headNE_.resize(A.num_row_, -1);
+
+  // temp stores the address of the latest entry in each row
+  std::vector<Int> temp(A.num_row_, -1);
+
+  // build linked lists for row access
+  for (Int col = 0; col < A.num_col_; ++col) {
+    for (Int el = A.start_[col]; el < A.start_[col + 1]; ++el) {
+      Int row = A.index_[el];
+      colNE_[el] = col;
+
+      if (temp[row] == -1)
+        headNE_[row] = el;
+      else
+        nextNE_[temp[row]] = el;
+
+      temp[row] = el;
+    }
+  }
+
   ptrNE_.clear();
   rowsNE_.clear();
 
@@ -164,34 +176,27 @@ Int FactorHiGHSSolver::buildNEstructureSparse(const HighsSparseMatrix& A,
 
     Int nz_in_col = 0;
 
-    for (Int rowEl = AT_.start_[row]; rowEl < AT_.start_[row + 1]; ++rowEl) {
-      Int col = AT_.index_[rowEl];
+    Int current = headNE_[row];
+    while (current != -1) {
+      Int col = colNE_[current];
 
-      // for each nonzero in the row, go down corresponding column
-      for (Int colEl = A.start_[col]; colEl < A.start_[col + 1]; ++colEl) {
+      // for each nonzero in the row, go down corresponding column, starting
+      // from current position
+      for (Int colEl = current; colEl < A.start_[col + 1]; ++colEl) {
         Int row2 = A.index_[colEl];
 
-        // skip when row2 is above row
-        if (row2 < row) continue;
+        // row2 is guaranteed to be larger or equal than row
+        // (provided that the columns of A are sorted)
 
         // save information that there is nonzero in position (row2,row).
         if (!is_nz[row2]) {
           is_nz[row2] = true;
           temp_index[nz_in_col] = row2;
           ++nz_in_col;
         }
-
-        // the same if statement could be executed without branching:
-        // (there does not seem to be a performance advantage)
-        /*
-          Int old_v = is_nz[row2];
-          Int new_v = 1;
-          Int diff_v = new_v - old_v;
-          is_nz[row2] = true;
-          temp_index[nz_in_col] = temp_index[nz_in_col] * old_v + diff_v * row2;
-          nz_in_col += diff_v;
-        */
       }
+
+      current = nextNE_[current];
     }
     // intersection of row with rows below finished.
 
@@ -213,89 +218,6 @@ Int FactorHiGHSSolver::buildNEstructureSparse(const HighsSparseMatrix& A,
   return kStatusOk;
 }
 
-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());
-
-  valNE_.assign(rowsNE_.size(), 0.0);
-
-  std::vector<double> work(A.num_row_, 0.0);
-
-  for (Int row = 0; row < A.num_row_; ++row) {
-    // go along the entries of the row, and then down each column.
-    // this builds the lower triangular part of the row-th column of AAt.
-    for (Int rowEl = AT_.start_[row]; rowEl < AT_.start_[row + 1]; ++rowEl) {
-      Int col = AT_.index_[rowEl];
-
-      const double theta =
-          scaling.empty() ? 1.0 : 1.0 / (scaling[col] + regul_.primal);
-
-      const double row_value = theta * AT_.value_[rowEl];
-
-      // for each nonzero in the row, go down corresponding column
-      for (Int colEl = A.start_[col]; colEl < A.start_[col + 1]; ++colEl) {
-        Int row2 = A.index_[colEl];
-
-        // skip when row2 is above row
-        if (row2 < row) continue;
-
-        // compute and accumulate value
-        double value = row_value * A.value_[colEl];
-        work[row2] += value;
-      }
-    }
-    // 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];
-      work[index] = 0.0;
-    }
-  }
-
-  return kStatusOk;
-}
-
-Int FactorHiGHSSolver::buildNEstructureJacek(const HighsSparseMatrix& A,
-                                             int64_t nz_limit) {
-  // Build lower triangular structure of AAt.
-  // This approach uses a column-wise copy of A and a collection of linked lists
-  // to access the rows
-
-  // NB: A must have sorted columns for this to work
-
-  // colNE stores the index of the column of each entry.
-  // It's useless when accessing the matrix by columns, but useful when doing it
-  // by rows.
-  colNE_.resize(A.numNz());
-
-  // head and next are a collection of linked lists to access the matrix by
-  // rows.
-  nextNE_.resize(A.numNz(), -1);
-  headNE_.resize(A.num_row_, -1);
-
-  // temp stores the address of the latest entry in each row
-  std::vector<Int> temp(A.num_row_, -1);
-
-  // build linked lists for row access
-  for (Int col = 0; col < A.num_col_; ++col) {
-    for (Int el = A.start_[col]; el < A.start_[col + 1]; ++el) {
-      Int row = A.index_[el];
-      colNE_[el] = col;
-
-      if (temp[row] == -1)
-        headNE_[row] = el;
-      else
-        nextNE_[temp[row]] = el;
-
-      temp[row] = el;
-    }
-  }
-}
-
 Int FactorHiGHSSolver::factorAS(const HighsSparseMatrix& A,
                                 const std::vector<double>& scaling) {
   // only execute factorisation if it has not been done yet
@@ -491,10 +413,8 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S, int64_t nz_limit) {
 
   log_.printDevInfo("Building NE structure\n");
 
-  buildNEstructureJacek(model_.A(), nz_limit);
-
   Clock clock;
-  if (Int status = buildNEstructureSparse(model_.A(), nz_limit)) return status;
+  if (Int status = buildNEstructure(model_.A(), nz_limit)) return status;
   if (info_) info_->matrix_structure_time = clock.stop();
 
   // create vector of signs of pivots
diff --git a/highs/ipm/hipo/ipm/FactorHiGHSSolver.h b/highs/ipm/hipo/ipm/FactorHiGHSSolver.h
@@ -37,16 +37,12 @@ class FactorHiGHSSolver : public LinearSolver {
   Int chooseNla();
   Int setNla();
   void setParallel();
-  Int buildNEstructureDense(const HighsSparseMatrix& A,
-                            int64_t nz_limit = kHighsIInf);
-  Int buildNEstructureSparse(const HighsSparseMatrix& A,
-                             int64_t nz_limit = kHighsIInf);
+
+  Int buildNEstructure(const HighsSparseMatrix& A,
+                       int64_t nz_limit = kHighsIInf);
   Int buildNEvalues(const HighsSparseMatrix& A,
                     const std::vector<double>& scaling);
 
-  Int buildNEstructureJacek(const HighsSparseMatrix& A,
-                            int64_t nz_limit = kHighsIInf);
-
   Int analyseAS(Symbolic& S);
   Int analyseNE(Symbolic& S, int64_t nz_limit = kHighsIInf);
 
