From: Denis Davydov <davydden@gmail.com>
Date: Tue, 9 Oct 2018 07:59:46 +0000 (+0200)
Subject: add DynamicSparsityPattern::compute_Tmmult_pattern() by moving existing code
X-Git-Tag: v9.1.0-rc1~638^2
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F7310%2Fhead;p=dealii.git

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 e188579eac..c33fe8edcb 100644
--- 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 2566bb0cb9..b7414ff150 100644
--- 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 37a1be647b..3caaa9d4a7 100644
--- 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 567340242a..62ca57ff73 100644
--- 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(