]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Fix quadratic behavior of SparsityPattern::copy_from().
authorWolfgang Bangerth <bangerth@colostate.edu>
Tue, 20 Mar 2018 23:06:31 +0000 (17:06 -0600)
committerWolfgang Bangerth <bangerth@colostate.edu>
Tue, 20 Mar 2018 23:06:31 +0000 (17:06 -0600)
source/lac/sparsity_pattern.cc

index 57663b5ff206db2252789a89f532baf4d352a453..4051dffa4c4d87f5dd3e94d73d1a14637769edf0 100644 (file)
@@ -545,13 +545,15 @@ void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
   // first init with the number of entries per row. if this matrix is square
   // then we also have to allocate memory for the diagonal entry, unless we
   // have already counted it
+  const bool matrix_is_square = (matrix.m() == matrix.n());
+
   std::vector<unsigned int> entries_per_row (matrix.m(), 0);
   for (size_type row=0; row<matrix.m(); ++row)
     {
       for (size_type col=0; col<matrix.n(); ++col)
         if (matrix(row,col) != 0)
           ++entries_per_row[row];
-      if ((matrix.m() == matrix.n())
+      if (matrix_is_square
           &&
           (matrix(row,row) == 0))
         ++entries_per_row[row];
@@ -559,11 +561,43 @@ void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
 
   reinit (matrix.m(), matrix.n(), entries_per_row);
 
-  // now set entries
+  // now set entries. if we enter entries row by row, then we'll get
+  // quadratic complexity in the number of entries per row. this is
+  // not usually a problem (we don't usually create dense matrices),
+  // but there are cases where it matters -- so we may as well be
+  // gentler and hand over a whole row of entries at a time
+  std::vector<size_type> column_indices;
+  column_indices.reserve (*std::max_element (entries_per_row.begin(),
+                                             entries_per_row.end()));
   for (size_type row=0; row<matrix.m(); ++row)
-    for (size_type col=0; col<matrix.n(); ++col)
-      if (matrix(row,col) != 0)
-        add (row,col);
+    {
+      column_indices.resize(entries_per_row[row]);
+
+      size_type current_index = 0;
+      for (size_type col=0; col<matrix.n(); ++col)
+        if (matrix(row,col) != 0)
+          {
+            column_indices[current_index] = col;
+            ++current_index;
+          }
+        else
+          // the (row,col) entry is zero; check if we need to add it
+          // anyway because it's the diagonal entry of a square
+          // matrix
+          if (matrix_is_square
+              &&
+              (col == row))
+            {
+              column_indices[current_index] = row;
+              ++current_index;
+            }
+
+      // check that we really added the correct number of indices
+      Assert (current_index == entries_per_row[row], ExcInternalError());
+
+      // now bulk add all of these entries
+      add_entries(row, column_indices.begin(), column_indices.end(), true);
+    }
 
   // finally compress
   compress ();

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.