NE structure now uses partial rw copy and corresponding entries

Author: filikat

highs/ipm/hipo/ipm/FactorHiGHSSolver.cpp

@@ -171,6 +171,8 @@ Int FactorHiGHSSolver::buildNEstructure(const HighsSparseMatrix& A,
   // temporary storage of indices
   std::vector<Int> temp_index(A.num_row_);
 
+  Clock clock;
+
   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.
@@ -216,6 +218,98 @@ Int FactorHiGHSSolver::buildNEstructure(const HighsSparseMatrix& A,
     }
   }
 
+  printf("time 1: %f\n", clock.stop());
+
+  return kStatusOk;
+}
+
+Int FactorHiGHSSolver::buildNEstructure_2(const HighsSparseMatrix& A,
+                                          int64_t nz_limit) {
+  // Build lower triangular structure of AAt.
+  // This approach uses a column-wise copy of A, a partial row-wise copy and a
+  // vector of corresponding indices.
+
+  // NB: A must have sorted columns for this to work
+
+  // create partial row-wise representation
+  // for now, using the full row-wise matrix
+  {
+    HighsSparseMatrix At = A;
+    At.ensureRowwise();
+    ptrNE_rw_ = std::move(At.start_);
+    idxNE_rw_ = std::move(At.index_);
+  }
+
+  {
+    std::vector<Int> temp = ptrNE_rw_;
+    corr_NE_.assign(A.numNz(), 0);
+
+    // build array or corresponding indices
+    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];
+
+        corr_NE_[temp[row]] = el;
+        temp[row]++;
+      }
+    }
+  }
+
+  std::vector<Int> ptrNE(A.num_row_ + 1);
+  std::vector<Int> rowsNE;
+
+  // keep track if given entry is nonzero, in column considered
+  std::vector<bool> is_nz(A.num_row_, false);
+
+  // temporary storage of indices
+  std::vector<Int> temp_index(A.num_row_);
+
+  Clock clock;
+
+  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.
+
+    Int nz_in_col = 0;
+
+    for (Int el = ptrNE_rw_[row]; el < ptrNE_rw_[row + 1]; ++el) {
+      Int col = idxNE_rw_[el];
+      Int corr = corr_NE_[el];
+
+      // for each nonzero in the row, go down corresponding column, starting
+      // from current position
+      for (Int colEl = corr; colEl < A.start_[col + 1]; ++colEl) {
+        Int row2 = A.index_[colEl];
+
+        // 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;
+        }
+      }
+    }
+    // intersection of row with rows below finished.
+
+    // if the total number of nonzeros exceeds the maximum, return error
+    if ((int64_t)ptrNE[row] + (int64_t)nz_in_col >= nz_limit) return kStatusOoM;
+
+    // update pointers
+    ptrNE[row + 1] = ptrNE[row] + nz_in_col;
+
+    // now assign indices
+    for (Int i = 0; i < nz_in_col; ++i) {
+      Int index = temp_index[i];
+      rowsNE.push_back(index);
+      is_nz[index] = false;
+    }
+  }
+
+  printf("time 2: %f\n", clock.stop());
+
   return kStatusOk;
 }
 
@@ -418,6 +512,8 @@ Int FactorHiGHSSolver::analyseNE(Symbolic& S, int64_t nz_limit) {
   if (Int status = buildNEstructure(model_.A(), nz_limit)) return status;
   if (info_) info_->matrix_structure_time = clock.stop();
 
+  buildNEstructure_2(model_.A(), nz_limit);
+
   // create vector of signs of pivots
   std::vector<Int> pivot_signs(model_.A().num_row_, 1);
 

highs/ipm/hipo/ipm/FactorHiGHSSolver.h

@@ -24,6 +24,9 @@ class FactorHiGHSSolver : public LinearSolver {
   std::vector<double> valNE_;
   std::vector<Int> colNE_, nextNE_, headNE_;
 
+  std::vector<Int> ptrNE_rw_, idxNE_rw_;
+  std::vector<Int> corr_NE_;
+
   const Regularisation& regul_;
 
   Info* info_ = nullptr;
@@ -39,6 +42,8 @@ class FactorHiGHSSolver : public LinearSolver {
 
   Int buildNEstructure(const HighsSparseMatrix& A,
                        int64_t nz_limit = kHighsIInf);
+  Int buildNEstructure_2(const HighsSparseMatrix& A,
+                         int64_t nz_limit = kHighsIInf);
   Int buildNEvalues(const HighsSparseMatrix& A,
                     const std::vector<double>& scaling);