Fix bad complexity in copy of sparsity pattern.

diff --git a/doc/news/changes.h b/doc/news/changes.h
index 690883db3040130b9f329180f9c091a1ad3cc1cb..cfef367ed01a68945fc2e6fe61326b979e75fa53 100644 (file)
--- a/doc/news/changes.h
+++ b/doc/news/changes.h
@@ -236,15 +236,15 @@ inconvenience this causes.
 
  <li>
  New: Added a new PolarManifold descriptor, that uses a polar coordinate
- system to compute new points, and modified the existing SphericalManifold 
- descriptor to use geodesics on the surface of the sphere. 
+ system to compute new points, and modified the existing SphericalManifold
+ descriptor to use geodesics on the surface of the sphere.
  <br>
  (Luca Heltai, Mauro Bardelloni, 2016/08/04)
  </li>
 
  <li>
- New: Added Python bindings to generate and manipulate a Triangulation from 
- Python. The Triangulation generated in Python can be saved and later, loaded 
+ New: Added Python bindings to generate and manipulate a Triangulation from
+ Python. The Triangulation generated in Python can be saved and later, loaded
  inside a C++ code.
  <br>
  (Bruno Turcksin, 2016/08/03)
@@ -495,6 +495,14 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+ <li> Improved: SparsityPattern::copy_from() copying from a
+ DynamicSparsityPattern argument had quadratic complexity in the number of
+ rows for sparsity patterns where most of the rows are of length zero. The bad
+ algorithm has been replaced by a linear complexity one.
+ <br>
+ (Dustin Kumor, Martin Kronbichler, 2016/08/31)
+ </li>
+
  <li> New: Rank-4 symmetric tensors of type SymmetricTensor can now
  be converted to rank-4 tensors of type Tensor.
  <br>

diff --git a/source/lac/sparsity_pattern.cc b/source/lac/sparsity_pattern.cc
index 501b4878ed507151427a204966ae60f920db498f..80fe41db7b36ea50a861c84f2c1401d990f36ea8 100644 (file)
--- a/source/lac/sparsity_pattern.cc
+++ b/source/lac/sparsity_pattern.cc
@@ -523,6 +523,46 @@ SparsityPattern::copy_from (const SparsityPatternType &dsp)
 
 
 
+// Use a special implementation for DynamicSparsityPattern where we can use
+// the column_number method to gain faster access to the
+// entries. DynamicSparsityPattern::iterator can show quadratic complexity in
+// case many rows are empty and the begin() method needs to jump to the next
+// free row. Otherwise, the code is exactly the same as above.
+template <>
+void
+SparsityPattern::copy_from (const DynamicSparsityPattern &dsp)
+{
+  const bool do_diag_optimize = (dsp.n_rows() == dsp.n_cols());
+  std::vector<unsigned int> row_lengths (dsp.n_rows());
+  for (size_type i=0; i<dsp.n_rows(); ++i)
+    {
+      row_lengths[i] = dsp.row_length(i);
+      if (do_diag_optimize && !dsp.exists(i,i))
+        ++row_lengths[i];
+    }
+  reinit (dsp.n_rows(), dsp.n_cols(), row_lengths);
+
+  if (n_rows() != 0 && n_cols() != 0)
+    for (size_type row = 0; row<dsp.n_rows(); ++row)
+      {
+        size_type *cols = &colnums[rowstart[row]] + (do_diag_optimize ? 1 : 0);
+        const unsigned int row_length = dsp.row_length(row);
+        for (unsigned int index=0; index < row_length; ++index)
+          {
+            const size_type col = dsp.column_number(row, index);
+            if ((col!=row) || !do_diag_optimize)
+              *cols++ = col;
+          }
+      }
+
+  // do not need to compress the sparsity pattern since we already have
+  // allocated the right amount of data, and the SparsityPatternType data is sorted,
+  // too.
+  compressed = true;
+}
+
+
+
 template <typename number>
 void SparsityPattern::copy_from (const FullMatrix<number> &matrix)
 {
@@ -998,7 +1038,6 @@ SparsityPattern::memory_consumption () const
 
 // explicit instantiations
 template void SparsityPattern::copy_from<SparsityPattern> (const SparsityPattern &);
-template void SparsityPattern::copy_from<DynamicSparsityPattern> (const DynamicSparsityPattern &);
 template void SparsityPattern::copy_from<float> (const FullMatrix<float> &);
 template void SparsityPattern::copy_from<double> (const FullMatrix<double> &);