From 2b61e4e9a7c0e25a2444db2d919f7e385e8e8e83 Mon Sep 17 00:00:00 2001
From: Christoph Goering <christoph.goering@tum.de>
Date: Fri, 15 Dec 2017 18:42:00 +0100
Subject: [PATCH] added DynamicSparsityPattern::compute_mmult_pattern(..)

---
 .../deal.II/lac/dynamic_sparsity_pattern.h    |  8 ++++
 source/lac/dynamic_sparsity_pattern.cc        | 45 +++++++++++++++++++
 2 files changed, 53 insertions(+)

diff --git a/include/deal.II/lac/dynamic_sparsity_pattern.h b/include/deal.II/lac/dynamic_sparsity_pattern.h
index 4d96e1d44a..b2ce716150 100644
--- a/include/deal.II/lac/dynamic_sparsity_pattern.h
+++ b/include/deal.II/lac/dynamic_sparsity_pattern.h
@@ -420,6 +420,14 @@ public:
    */
   void symmetrize ();
 
+  /**
+   * Construct and store in this object the sparsity pattern corresponding to
+   * the product of @p left and @p right sparsity pattern.
+   */
+  template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
+  void compute_mmult_pattern(const SparsityPatternTypeLeft &left,
+                             const SparsityPatternTypeRight &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/source/lac/dynamic_sparsity_pattern.cc b/source/lac/dynamic_sparsity_pattern.cc
index 18fce9577e..e41f1c05de 100644
--- a/source/lac/dynamic_sparsity_pattern.cc
+++ b/source/lac/dynamic_sparsity_pattern.cc
@@ -14,6 +14,7 @@
 // ---------------------------------------------------------------------
 
 #include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/sparsity_pattern.h>
 #include <deal.II/base/memory_consumption.h>
 
 #include <algorithm>
@@ -407,6 +408,37 @@ DynamicSparsityPattern::symmetrize ()
 
 
 
+template <typename SparsityPatternTypeLeft, typename SparsityPatternTypeRight>
+void
+DynamicSparsityPattern::compute_mmult_pattern (
+  const SparsityPatternTypeLeft &left,
+  const SparsityPatternTypeRight &right)
+{
+  Assert (left.n_cols() == right.n_rows(),
+          ExcDimensionMismatch(left.n_cols(), right.n_rows()));
+
+  this->reinit(left.n_rows(),right.n_cols());
+
+  typename SparsityPatternTypeLeft::iterator
+  it_left = left.begin(),
+  end_left = left.end();
+  for (; it_left!=end_left; ++it_left)
+    {
+      const unsigned int j = it_left->column();
+
+      // We are sitting on entry (i,j) of the left sparsity pattern. We then need
+      // to add all entries (i,k) to the final sparsity pattern where (j,k) exists
+      // in the right sparsity pattern -- i.e., we need to iterate over row j.
+      typename SparsityPatternTypeRight::iterator
+      it_right = right.begin(j),
+      end_right = right.end(j);
+      for (; it_right!=end_right; ++it_right)
+        this->add(it_left->row(),it_right->column());
+    }
+}
+
+
+
 void
 DynamicSparsityPattern::print (std::ostream &out) const
 {
@@ -519,4 +551,17 @@ add_entries(std::vector<size_type>::iterator,
             const bool);
 #endif
 
+template void DynamicSparsityPattern::compute_mmult_pattern(
+  const DynamicSparsityPattern &,
+  const DynamicSparsityPattern &);
+template void DynamicSparsityPattern::compute_mmult_pattern(
+  const DynamicSparsityPattern &,
+  const SparsityPattern &);
+template void DynamicSparsityPattern::compute_mmult_pattern(
+  const SparsityPattern &,
+  const DynamicSparsityPattern &);
+template void DynamicSparsityPattern::compute_mmult_pattern(
+  const SparsityPattern &,
+  const SparsityPattern &);
+
 DEAL_II_NAMESPACE_CLOSE
-- 
2.39.5