apply requested adjustments

resolve/matrix/Utilities.cpp

@@ -79,16 +79,16 @@ namespace ReSolve
       std::unique_ptr<index_type[]> used(new index_type[n_rows]);
       std::fill_n(used.get(), n_rows, 0);
 
-      // allocation stage
-
-      // NOTE: so aside from tracking the number of mirrored values for allocation
-      //       purposes, we also have to compute the amount of nonzeroes in each row
-      //       since we're dealing with a COO input matrix. we store these values in
-      //       row + 1 slots in new_rows, to avoid allocating extra space :)
+      // allocation stage, the first loop is O(nnz) no matter the branch and the second
+      // is O(n)
+      //
+      // all this does is prepare the row index array based on the nonzeroes stored in
+      // the input coo matrix. if the matrix is symmetric and either upper or lower
+      // triangular, it additionally corrects for mirrored values using the `used` array
+      // and validates the triangularity. the branch is done to avoid the extra work if
+      // it's not necessary
 
       if (!A->symmetric() || A->expanded()) {
-        // NOTE: this branch is special cased because there is a bunch of extra
-        //       bookkeeping done if a matrix is symmetric and unexpanded
         for (index_type i = 0; i < nnz; i++) {
           new_rows[rows[i] + 1]++;
         }
@@ -114,23 +114,27 @@ namespace ReSolve
       }
 
       index_type* new_columns = new index_type[new_rows[n_rows]];
-      // TODO: we need to find a better way to fix this. what this fixes is
-      //       that otherwise, the column values are all zero by default and
-      //       this seems to mess with the way the insertion position is
-      //       selected if the column of the pair we're looking to insert
-      //       is zero
       std::fill_n(new_columns, new_rows[n_rows], -1);
       real_type* new_values = new real_type[new_rows[n_rows]];
 
-      // fill stage
+      // fill stage, approximately O(nnz * m) in the worst case
+      //
+      // all this does is iterate over the nonzeroes in the coo matrix,
+      // check to see if a value at that colum already exists using binary search,
+      // and if it does, then insert the new value at that position (deduplicating
+      // the matrix), otherwise, it allocates a new spot in the row (where you see
+      // used[rows[i]]++) and shifts everything over, performing what is effectively
+      // insertion sort. the lower half is conditioned on the matrix being symmetric
+      // and stored as either upper-triangular or lower-triangular, and just
+      // performs the same as what is described above, but with the indices swapped.
 
       for (index_type i = 0; i < nnz; i++) {
         index_type insertion_pos =
             static_cast<index_type>(
                 std::lower_bound(&new_columns[new_rows[rows[i]]],
                                  &new_columns[new_rows[rows[i]] + used[rows[i]]],
-                                 columns[i]) -
-                new_columns);
+                                 columns[i])
+                - new_columns);
 
         // NOTE: the stuff beyond here is basically insertion sort. i'm not
         //       sure if it's the best way of going about sorting the
@@ -163,8 +167,8 @@ namespace ReSolve
               static_cast<index_type>(
                   std::lower_bound(&new_columns[new_rows[columns[i]]],
                                    &new_columns[new_rows[columns[i]] + used[columns[i]]],
-                                   rows[i]) -
-                  new_columns);
+                                   rows[i])
+                  - new_columns);
 
           if (new_columns[mirrored_insertion_pos] == rows[i]) {
             new_values[mirrored_insertion_pos] = values[i];
@@ -182,7 +186,13 @@ namespace ReSolve
         }
       }
 
-      // backshifting stage
+      // backshifting stage, approximately O(n^2 * nnz) in the worst case
+      // (this probably isn't the tightest bound i could apply)
+      //
+      // all this does is shift back rows to remove empty space in between them
+      // by indexing each row in order, checking to see if the amount of used
+      // spaces is equivalent to the amount of nonzeroes in the row, and if not,
+      // shifts every subsequent row back the amount of unused spaces
 
       for (index_type row = 0; row < n_rows - 1; row++) {
         index_type row_nnz = new_rows[row + 1] - new_rows[row];
@@ -206,18 +216,25 @@ namespace ReSolve
         }
       }
 
-      if (B->destroyMatrixData(memory::HOST) != 0 ||
-          B->setMatrixData(new_rows, new_columns, new_values, memory::HOST) != 0) {
+      B->setSymmetric(A->symmetric());
+      B->setNnz(new_rows[n_rows]);
+      // NOTE: this is necessary because updateData always reads the current nnz from
+      //       this field
+      B->setNnzExpanded(new_rows[n_rows]);
+      B->setExpanded(true);
+
+      if (B->updateData(new_rows, new_columns, new_values, memory::HOST, memspace) != 0) {
+        delete[] new_rows;
+        delete[] new_columns;
+        delete[] new_values;
+
         assert(false && "invalid state after coo -> csr conversion");
         return -1;
       }
 
-      // TODO: set data ownership / value ownership here. we'd be leaking
-      //       memory here otherwise
-
-      B->setSymmetric(A->symmetric());
-      B->setNnz(new_rows[n_rows]);
-      B->setExpanded(true);
+      delete[] new_rows;
+      delete[] new_columns;
+      delete[] new_values;
 
       return 0;
     }