Improve inverse connectivity construction efficiency

One does not use std::map anymore

src/mesh/ConnectivityComputer.cpp

@@ -4,8 +4,6 @@
 #include <utils/Exceptions.hpp>
 
 #include <iostream>
-#include <map>
-#include <sstream>
 
 template <typename ConnectivityType>
 PUGS_INLINE ConnectivityMatrix
@@ -17,8 +15,6 @@ ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connecti
   if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
     const ConnectivityMatrix& child_to_item_matrix = connectivity.getMatrix(child_item_type, item_type);
 
-    std::cout << "computing connectivity " << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
-
     item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix);
   } else {
     std::stringstream error_msg;
@@ -44,32 +40,46 @@ template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(cons
 ConnectivityMatrix
 ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
 {
-  std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map;
-  const size_t& number_of_items = item_to_child_matrix.numberOfRows();
+  const size_t& number_of_rows = item_to_child_matrix.numberOfRows();
 
-  for (unsigned int j = 0; j < number_of_items; ++j) {
-    const auto& item_to_child = item_to_child_matrix[j];
-    for (unsigned int r = 0; r < item_to_child.size(); ++r) {
-      child_to_item_vector_map[item_to_child[r]].push_back(j);
-    }
-  }
+  if ((item_to_child_matrix.values().size() > 0)) {
+    const size_t& number_of_columns = max(item_to_child_matrix.values());
 
-  {
-    size_t i = 0;
-    for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
-      if (child_item_id != i) {
-        throw NotImplementedError("sparse item numerotation");
+    Array<uint32_t> transposed_next_free_column_index(number_of_columns + 1);
+    transposed_next_free_column_index.fill(0);
+
+    Array<uint32_t> transposed_rows_map(transposed_next_free_column_index.size() + 1);
+    transposed_rows_map.fill(0);
+    for (size_t i = 0; i < number_of_rows; ++i) {
+      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
+        transposed_rows_map[item_to_child_matrix.values()[j] + 1]++;
       }
-      ++i;
     }
+    for (size_t i = 1; i < transposed_rows_map.size(); ++i) {
+      transposed_rows_map[i] += transposed_rows_map[i - 1];
     }
+    Array<uint32_t> transposed_column_indices(transposed_rows_map[transposed_rows_map.size() - 1]);
 
-  std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size());
-  for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
-    child_to_items_vector[child_item_id] = item_vector;
+    for (size_t i = 0; i < number_of_rows; ++i) {
+      for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
+        size_t i_column_index = item_to_child_matrix.values()[j];
+        uint32_t& shift       = transposed_next_free_column_index[i_column_index];
+
+        transposed_column_indices[transposed_rows_map[i_column_index] + shift] = i;
+        Assert(shift < transposed_rows_map[i_column_index + 1]);
+        shift++;
+      }
     }
 
-  return ConnectivityMatrix(child_to_items_vector);
+    return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
+  } else {
+    // empty connectivity
+    Array<uint32_t> transposed_rows_map{1};
+    transposed_rows_map[0] = 0;
+    Array<uint32_t> transposed_column_indices;
+
+    return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
+  }
 }
 
 template <typename ItemOfItem, typename ConnectivityType>

src/mesh/ConnectivityMatrix.hpp

@@ -56,7 +56,22 @@ class ConnectivityMatrix
   }
 
   PUGS_INLINE
-  ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) : m_is_built{true}
+  ConnectivityMatrix(const Array<const uint32_t>& row_map,
+                     const Array<const uint32_t>& column_indices) noexcept(NO_ASSERT)
+    : m_row_map{row_map}, m_column_indices{column_indices}, m_is_built{true}
+  {
+    Assert(m_column_indices.size() == m_row_map[m_row_map.size() - 1], "incompatible row map and column indices");
+    Assert(m_row_map.size() > 0, "invalid row map");
+    Assert(m_row_map[0] == 0, "row map should start with 0");
+#ifndef NDEBUG
+    for (size_t i = 1; i < m_row_map.size(); ++i) {
+      Assert(m_row_map[i] > m_row_map[i - 1], "row map values must be strictly increasing");
+    }
+#endif   // NDEBUG
+  }
+
+  PUGS_INLINE
+  ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) noexcept : m_is_built{true}
   {
     m_row_map = [&] {
       Array<uint32_t> row_map(initializer.size() + 1);