add DynamicSparsityPattern::compute_Tmmult_pattern() by moving existing code

diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h
index e188579eac3a35cf0323bc56b3601959f5549e6a..c33fe8edcbe34b33e76a1566cd727a8d86ef3543 100644 (file)
--- a/include/deal.II/lac/dynamic_sparsity_pattern.h
+++ b/include/deal.II/lac/dynamic_sparsity_pattern.h
@@ -32,7 +32,7 @@
 DEAL_II_NAMESPACE_OPEN
 
 class DynamicSparsityPattern;
-
+class SparsityPattern;
 
 /*! @addtogroup Sparsity
  *@{
@@ -446,6 +446,14 @@ public:
   compute_mmult_pattern(const SparsityPatternTypeLeft & left,
                         const SparsityPatternTypeRight &right);
 
+  /**
+   * Construct and store in this object the sparsity pattern corresponding to
+   * the product of transposed @p left and and non-transpose @p right sparsity pattern.
+   */
+  void
+  compute_Tmmult_pattern(const SparsityPattern &left,
+                         const SparsityPattern &right);
+
   /**
    * Print the sparsity pattern. The output consists of one line per row of
    * the format <tt>[i,j1,j2,j3,...]</tt>. <i>i</i> is the row number and

diff --git a/include/deal.II/lac/sparse_matrix.templates.h b/include/deal.II/lac/sparse_matrix.templates.h
index 2566bb0cb9a60a7ac45b3e4280d3e132e72c252f..b7414ff15025b86ba2b455ffcb4035dd70c24171 100644 (file)
--- a/include/deal.II/lac/sparse_matrix.templates.h
+++ b/include/deal.II/lac/sparse_matrix.templates.h
@@ -1145,39 +1145,10 @@ SparseMatrix<number>::Tmmult(SparseMatrix<numberC> &      C,
       C.clear();
       sp_C.reinit(0, 0, 0);
 
-      // create a sparsity pattern for the matrix. we will go through all the
-      // rows in the matrix A, and for each column in a row we add the whole
-      // row of matrix B with that row number. This means that we will insert
-      // a lot of entries to each row, which is best handled by the
-      // DynamicSparsityPattern class.
+      // create a sparsity pattern for the matrix.
       {
-        DynamicSparsityPattern dsp(n(), B.n());
-        for (size_type i = 0; i < sp_A.n_rows(); ++i)
-          {
-            const size_type *      rows = &sp_A.colnums[sp_A.rowstart[i]];
-            const size_type *const end_rows =
-              &sp_A.colnums[sp_A.rowstart[i + 1]];
-            // cast away constness to conform with dsp.add_entries interface
-            size_type *new_cols =
-              const_cast<size_type *>(&sp_B.colnums[sp_B.rowstart[i]]);
-            size_type *end_new_cols =
-              const_cast<size_type *>(&sp_B.colnums[sp_B.rowstart[i + 1]]);
-
-            if (sp_B.n_rows() == sp_B.n_cols())
-              ++new_cols;
-
-            for (; rows != end_rows; ++rows)
-              {
-                const size_type row = *rows;
-
-                // if B has a diagonal, need to add that manually. this way,
-                // we maintain sortedness.
-                if (sp_B.n_rows() == sp_B.n_cols())
-                  dsp.add(row, i);
-
-                dsp.add_entries(row, new_cols, end_new_cols, true);
-              }
-          }
+        DynamicSparsityPattern dsp;
+        dsp.compute_Tmmult_pattern(sp_A, sp_B);
         sp_C.copy_from(dsp);
       }
 

diff --git a/include/deal.II/lac/sparsity_pattern.h b/include/deal.II/lac/sparsity_pattern.h
index 37a1be647bb282d413b2526b126b3d07219c64bd..3caaa9d4a78e71565d7d308958668b419153fb44 100644 (file)
--- a/include/deal.II/lac/sparsity_pattern.h
+++ b/include/deal.II/lac/sparsity_pattern.h
@@ -1129,6 +1129,7 @@ private:
   friend class ChunkSparseMatrix;
 
   friend class ChunkSparsityPattern;
+  friend class DynamicSparsityPattern;
 
   /**
    * Also give access to internal details to the iterator/accessor classes.

diff --git a/source/lac/dynamic_sparsity_pattern.cc b/source/lac/dynamic_sparsity_pattern.cc
index 567340242a6ca1b0300010317f4a77c36f2ba980..62ca57ff7375616c47dc97e0d2af68c70848ba43 100644 (file)
--- a/source/lac/dynamic_sparsity_pattern.cc
+++ b/source/lac/dynamic_sparsity_pattern.cc
@@ -405,6 +405,48 @@ DynamicSparsityPattern::symmetrize()
 
 
 
+void
+DynamicSparsityPattern::compute_Tmmult_pattern(const SparsityPattern &sp_A,
+                                               const SparsityPattern &sp_B)
+{
+  Assert(sp_A.n_rows() == sp_B.n_rows(),
+         ExcDimensionMismatch(sp_A.n_rows(), sp_B.n_rows()));
+
+  this->reinit(sp_A.n_cols(), sp_B.n_cols());
+  // we will go through all the
+  // rows in the matrix A, and for each column in a row we add the whole
+  // row of matrix B with that row number. This means that we will insert
+  // a lot of entries to each row, which is best handled by the
+  // DynamicSparsityPattern class.
+  for (size_type i = 0; i < sp_A.n_rows(); ++i)
+    {
+      const size_type *      rows     = &sp_A.colnums[sp_A.rowstart[i]];
+      const size_type *const end_rows = &sp_A.colnums[sp_A.rowstart[i + 1]];
+      // cast away constness to conform with dsp.add_entries interface
+      size_type *new_cols =
+        const_cast<size_type *>(&sp_B.colnums[sp_B.rowstart[i]]);
+      size_type *end_new_cols =
+        const_cast<size_type *>(&sp_B.colnums[sp_B.rowstart[i + 1]]);
+
+      if (sp_B.n_rows() == sp_B.n_cols())
+        ++new_cols;
+
+      for (; rows != end_rows; ++rows)
+        {
+          const size_type row = *rows;
+
+          // if B has a diagonal, need to add that manually. this way,
+          // we maintain sortedness.
+          if (sp_B.n_rows() == sp_B.n_cols())
+            this->add(row, i);
+
+          this->add_entries(row, new_cols, end_new_cols, true);
+        }
+    }
+}
+
+
+
 template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
 void
 DynamicSparsityPattern::compute_mmult_pattern(