Precompile TensorProductMatrixSymmetricSum kernels

diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h
index ffaa3de04c97dc34dbb749d9a106024ce05869a2..0d2ef2a199d55edf652f86bab12e41135f7cf109 100644 (file)
--- a/include/deal.II/lac/tensor_product_matrix.h
+++ b/include/deal.II/lac/tensor_product_matrix.h
@@ -21,6 +21,7 @@
 
 #include <deal.II/base/array_view.h>
 #include <deal.II/base/mutex.h>
+#include <deal.II/base/vectorization.h>
 
 #include <deal.II/lac/lapack_full_matrix.h>
 
@@ -502,6 +503,167 @@ TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::n() const
 }
 
 
+namespace internal
+{
+  namespace TensorProductMatrixSymmetricSum
+  {
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    vmult(Number *                                 dst,
+          const Number *                           src,
+          AlignedVector<Number> &                  tmp,
+          const std::array<Table<2, Number>, dim> &mass_matrix,
+          const std::array<Table<2, Number>, dim> &derivative_matrix)
+    {
+      const unsigned int n_rows_1d = mass_matrix[0].n_rows();
+      const unsigned int n         = Utilities::fixed_power<dim>(n_rows_1d);
+
+      tmp.resize_fast(n * 2);
+      Number *t = tmp.begin();
+
+      internal::EvaluatorTensorProduct<internal::evaluate_general,
+                                       dim,
+                                       n_rows_1d_templated,
+                                       n_rows_1d_templated,
+                                       Number>
+        eval({}, {}, {}, n_rows_1d, n_rows_1d);
+
+      if (dim == 1)
+        {
+          const Number *A = &derivative_matrix[0](0, 0);
+          eval.template apply<0, false, false>(A, src, dst);
+        }
+
+      else if (dim == 2)
+        {
+          const Number *A0 = &derivative_matrix[0](0, 0);
+          const Number *M0 = &mass_matrix[0](0, 0);
+          const Number *A1 = &derivative_matrix[1](0, 0);
+          const Number *M1 = &mass_matrix[1](0, 0);
+          eval.template apply<0, false, false>(M0, src, t);
+          eval.template apply<1, false, false>(A1, t, dst);
+          eval.template apply<0, false, false>(A0, src, t);
+          eval.template apply<1, false, true>(M1, t, dst);
+        }
+
+      else if (dim == 3)
+        {
+          const Number *A0 = &derivative_matrix[0](0, 0);
+          const Number *M0 = &mass_matrix[0](0, 0);
+          const Number *A1 = &derivative_matrix[1](0, 0);
+          const Number *M1 = &mass_matrix[1](0, 0);
+          const Number *A2 = &derivative_matrix[2](0, 0);
+          const Number *M2 = &mass_matrix[2](0, 0);
+          eval.template apply<0, false, false>(M0, src, t + n);
+          eval.template apply<1, false, false>(M1, t + n, t);
+          eval.template apply<2, false, false>(A2, t, dst);
+          eval.template apply<1, false, false>(A1, t + n, t);
+          eval.template apply<0, false, false>(A0, src, t + n);
+          eval.template apply<1, false, true>(M1, t + n, t);
+          eval.template apply<2, false, true>(M2, t, dst);
+        }
+
+      else
+        AssertThrow(false, ExcNotImplemented());
+    }
+
+
+
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    apply_inverse(Number *                                      dst,
+                  const Number *                                src,
+                  AlignedVector<Number> &                       tmp,
+                  const std::array<Table<2, Number>, dim> &     eigenvectors,
+                  const std::array<AlignedVector<Number>, dim> &eigenvalues)
+    {
+      const unsigned int n_rows_1d = eigenvectors[0].n_rows();
+      const unsigned int n         = Utilities::fixed_power<dim>(n_rows_1d);
+
+      tmp.resize_fast(n);
+      Number *t = tmp.begin();
+
+      internal::EvaluatorTensorProduct<internal::evaluate_general,
+                                       dim,
+                                       n_rows_1d_templated,
+                                       n_rows_1d_templated,
+                                       Number>
+        eval({}, {}, {}, n_rows_1d, n_rows_1d);
+
+      // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
+      //       --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
+      //       --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
+      //       while the eigenvectors are stored column-wise in S, i.e.
+      //       rows correspond to dofs whereas columns to eigenvalue indices!
+      if (dim == 1)
+        {
+          const Number *S = &eigenvectors[0](0, 0);
+          eval.template apply<0, true, false>(S, src, t);
+          for (unsigned int i = 0; i < n_rows_1d; ++i)
+            t[i] /= eigenvalues[0][i];
+          eval.template apply<0, false, false>(S, t, dst);
+        }
+
+      else if (dim == 2)
+        {
+          const Number *S0 = &(eigenvectors[0](0, 0));
+          const Number *S1 = &(eigenvectors[1](0, 0));
+          eval.template apply<0, true, false>(S0, src, t);
+          eval.template apply<1, true, false>(S1, t, dst);
+          for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1)
+            for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
+              dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
+          eval.template apply<0, false, false>(S0, dst, t);
+          eval.template apply<1, false, false>(S1, t, dst);
+        }
+
+      else if (dim == 3)
+        {
+          const Number *S0 = &eigenvectors[0](0, 0);
+          const Number *S1 = &eigenvectors[1](0, 0);
+          const Number *S2 = &eigenvectors[2](0, 0);
+          eval.template apply<0, true, false>(S0, src, t);
+          eval.template apply<1, true, false>(S1, t, dst);
+          eval.template apply<2, true, false>(S2, dst, t);
+          for (unsigned int i2 = 0, c = 0; i2 < n_rows_1d; ++i2)
+            for (unsigned int i1 = 0; i1 < n_rows_1d; ++i1)
+              for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c)
+                t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] +
+                         eigenvalues[0][i0]);
+          eval.template apply<0, false, false>(S0, t, dst);
+          eval.template apply<1, false, false>(S1, dst, t);
+          eval.template apply<2, false, false>(S2, t, dst);
+        }
+
+      else
+        Assert(false, ExcNotImplemented());
+    }
+
+
+
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    select_vmult(Number *                                 dst,
+                 const Number *                           src,
+                 AlignedVector<Number> &                  tmp,
+                 const std::array<Table<2, Number>, dim> &mass_matrix,
+                 const std::array<Table<2, Number>, dim> &derivative_matrix);
+
+
+
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    select_apply_inverse(
+      Number *                                      dst,
+      const Number *                                src,
+      AlignedVector<Number> &                       tmp,
+      const std::array<Table<2, Number>, dim> &     eigenvectors,
+      const std::array<AlignedVector<Number>, dim> &eigenvalues);
+  } // namespace TensorProductMatrixSymmetricSum
+
+} // namespace internal
+
+
 
 template <int dim, typename Number, int n_rows_1d>
 inline void
@@ -512,61 +674,17 @@ TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::vmult(
   AssertDimension(dst_view.size(), this->m());
   AssertDimension(src_view.size(), this->n());
   std::lock_guard<std::mutex> lock(this->mutex);
-  const unsigned int          n = Utilities::fixed_power<dim>(
-    n_rows_1d > 0 ? n_rows_1d : eigenvalues[0].size());
-  tmp_array.resize_fast(n * 2);
-  constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0;
-  internal::EvaluatorTensorProduct<internal::evaluate_general,
-                                   dim,
-                                   kernel_size,
-                                   kernel_size,
-                                   Number>
-                eval(AlignedVector<Number>{},
-         AlignedVector<Number>{},
-         AlignedVector<Number>{},
-         mass_matrix[0].n_rows(),
-         mass_matrix[0].n_rows());
-  Number *      t   = tmp_array.begin();
-  const Number *src = src_view.begin();
-  Number *      dst = dst_view.data();
 
-  if (dim == 1)
-    {
-      const Number *A = &derivative_matrix[0](0, 0);
-      eval.template apply<0, false, false>(A, src, dst);
-    }
-
-  else if (dim == 2)
-    {
-      const Number *A0 = &derivative_matrix[0](0, 0);
-      const Number *M0 = &mass_matrix[0](0, 0);
-      const Number *A1 = &derivative_matrix[1](0, 0);
-      const Number *M1 = &mass_matrix[1](0, 0);
-      eval.template apply<0, false, false>(M0, src, t);
-      eval.template apply<1, false, false>(A1, t, dst);
-      eval.template apply<0, false, false>(A0, src, t);
-      eval.template apply<1, false, true>(M1, t, dst);
-    }
-
-  else if (dim == 3)
-    {
-      const Number *A0 = &derivative_matrix[0](0, 0);
-      const Number *M0 = &mass_matrix[0](0, 0);
-      const Number *A1 = &derivative_matrix[1](0, 0);
-      const Number *M1 = &mass_matrix[1](0, 0);
-      const Number *A2 = &derivative_matrix[2](0, 0);
-      const Number *M2 = &mass_matrix[2](0, 0);
-      eval.template apply<0, false, false>(M0, src, t + n);
-      eval.template apply<1, false, false>(M1, t + n, t);
-      eval.template apply<2, false, false>(A2, t, dst);
-      eval.template apply<1, false, false>(A1, t + n, t);
-      eval.template apply<0, false, false>(A0, src, t + n);
-      eval.template apply<1, false, true>(M1, t + n, t);
-      eval.template apply<2, false, true>(M2, t, dst);
-    }
+  Number *      dst = dst_view.begin();
+  const Number *src = src_view.begin();
 
+  if (n_rows_1d != -1)
+    internal::TensorProductMatrixSymmetricSum::vmult<
+      n_rows_1d == -1 ? 0 : n_rows_1d>(
+      dst, src, tmp_array, mass_matrix, derivative_matrix);
   else
-    AssertThrow(false, ExcNotImplemented());
+    internal::TensorProductMatrixSymmetricSum::select_vmult<1>(
+      dst, src, tmp_array, mass_matrix, derivative_matrix);
 }
 
 
@@ -592,70 +710,17 @@ TensorProductMatrixSymmetricSumBase<dim, Number, n_rows_1d>::apply_inverse(
 {
   AssertDimension(dst_view.size(), this->n());
   AssertDimension(src_view.size(), this->m());
-  const unsigned int n = n_rows_1d > 0 ? n_rows_1d : eigenvalues[0].size();
-  tmp_array.resize_fast(Utilities::fixed_power<dim>(n));
-  constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0;
-  internal::EvaluatorTensorProduct<internal::evaluate_general,
-                                   dim,
-                                   kernel_size,
-                                   kernel_size,
-                                   Number>
-                eval(AlignedVector<Number>(),
-         AlignedVector<Number>(),
-         AlignedVector<Number>(),
-         mass_matrix[0].n_rows(),
-         mass_matrix[0].n_rows());
-  Number *      t   = tmp_array.begin();
-  const Number *src = src_view.data();
-  Number *      dst = dst_view.data();
-
-  // NOTE: dof_to_quad has to be interpreted as 'dof to eigenvalue index'
-  //       --> apply<.,true,.> (S,src,dst) calculates dst = S^T * src,
-  //       --> apply<.,false,.> (S,src,dst) calculates dst = S * src,
-  //       while the eigenvectors are stored column-wise in S, i.e.
-  //       rows correspond to dofs whereas columns to eigenvalue indices!
-  if (dim == 1)
-    {
-      const Number *S = &eigenvectors[0](0, 0);
-      eval.template apply<0, true, false>(S, src, t);
-      for (unsigned int i = 0; i < n; ++i)
-        t[i] /= eigenvalues[0][i];
-      eval.template apply<0, false, false>(S, t, dst);
-    }
-
-  else if (dim == 2)
-    {
-      const Number *S0 = &(eigenvectors[0](0, 0));
-      const Number *S1 = &(eigenvectors[1](0, 0));
-      eval.template apply<0, true, false>(S0, src, t);
-      eval.template apply<1, true, false>(S1, t, dst);
-      for (unsigned int i1 = 0, c = 0; i1 < n; ++i1)
-        for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
-          dst[c] /= (eigenvalues[1][i1] + eigenvalues[0][i0]);
-      eval.template apply<0, false, false>(S0, dst, t);
-      eval.template apply<1, false, false>(S1, t, dst);
-    }
 
-  else if (dim == 3)
-    {
-      const Number *S0 = &eigenvectors[0](0, 0);
-      const Number *S1 = &eigenvectors[1](0, 0);
-      const Number *S2 = &eigenvectors[2](0, 0);
-      eval.template apply<0, true, false>(S0, src, t);
-      eval.template apply<1, true, false>(S1, t, dst);
-      eval.template apply<2, true, false>(S2, dst, t);
-      for (unsigned int i2 = 0, c = 0; i2 < n; ++i2)
-        for (unsigned int i1 = 0; i1 < n; ++i1)
-          for (unsigned int i0 = 0; i0 < n; ++i0, ++c)
-            t[c] /=
-              (eigenvalues[2][i2] + eigenvalues[1][i1] + eigenvalues[0][i0]);
-      eval.template apply<0, false, false>(S0, t, dst);
-      eval.template apply<1, false, false>(S1, dst, t);
-      eval.template apply<2, false, false>(S2, t, dst);
-    }
+  Number *      dst = dst_view.begin();
+  const Number *src = src_view.begin();
 
+  if (n_rows_1d != -1)
+    internal::TensorProductMatrixSymmetricSum::apply_inverse<
+      n_rows_1d == -1 ? 0 : n_rows_1d>(
+      dst, src, tmp_array, eigenvectors, eigenvalues);
   else
-    Assert(false, ExcNotImplemented());
+    internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>(
+      dst, src, tmp_array, eigenvectors, eigenvalues);
 }
 
 

diff --git a/include/deal.II/lac/tensor_product_matrix.templates.h b/include/deal.II/lac/tensor_product_matrix.templates.h
new file mode 100644 (file)
index 0000000..4bbb26f
--- /dev/null
+++ b/include/deal.II/lac/tensor_product_matrix.templates.h
@@ -0,0 +1,86 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#ifndef dealii_tensor_product_matrix_templates_h
+#define dealii_tensor_product_matrix_templates_h
+
+
+#include <deal.II/base/config.h>
+
+#include <deal.II/lac/tensor_product_matrix.h>
+
+#ifndef FDM_DEGREE_MAX
+// We set this value 17, since this is the value needed for smoothers for
+// cell-centered patches with overlap for continuous elements with degrees up
+// to FE_EVAL_FACTORY_DEGREE_MAX.
+#  define FDM_DEGREE_MAX 17
+#endif
+
+DEAL_II_NAMESPACE_OPEN
+
+
+namespace internal
+{
+  namespace TensorProductMatrixSymmetricSum
+  {
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    select_vmult(Number *                                 dst,
+                 const Number *                           src,
+                 AlignedVector<Number> &                  tmp,
+                 const std::array<Table<2, Number>, dim> &mass_matrix,
+                 const std::array<Table<2, Number>, dim> &derivative_matrix)
+    {
+      const int n_rows_1d = mass_matrix[0].n_rows();
+
+      if (n_rows_1d_templated == n_rows_1d)
+        vmult<n_rows_1d_templated>(
+          dst, src, tmp, mass_matrix, derivative_matrix);
+      else if (n_rows_1d_templated < FDM_DEGREE_MAX)
+        select_vmult<std::min(n_rows_1d_templated + 1, FDM_DEGREE_MAX)>(
+          dst, src, tmp, mass_matrix, derivative_matrix);
+      else
+        vmult<0>(dst, src, tmp, mass_matrix, derivative_matrix);
+    }
+
+
+
+    template <int n_rows_1d_templated, std::size_t dim, typename Number>
+    void
+    select_apply_inverse(
+      Number *                                      dst,
+      const Number *                                src,
+      AlignedVector<Number> &                       tmp,
+      const std::array<Table<2, Number>, dim> &     eigenvectors,
+      const std::array<AlignedVector<Number>, dim> &eigenvalues)
+    {
+      const int n_rows_1d = eigenvectors[0].n_rows();
+
+      if (n_rows_1d_templated == n_rows_1d)
+        apply_inverse<n_rows_1d_templated>(
+          dst, src, tmp, eigenvectors, eigenvalues);
+      else if (n_rows_1d_templated < FDM_DEGREE_MAX)
+        select_apply_inverse<std::min(n_rows_1d_templated + 1, FDM_DEGREE_MAX)>(
+          dst, src, tmp, eigenvectors, eigenvalues);
+      else
+        apply_inverse<0>(dst, src, tmp, eigenvectors, eigenvalues);
+    }
+  } // namespace TensorProductMatrixSymmetricSum
+} // namespace internal
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif

diff --git a/source/lac/CMakeLists.txt b/source/lac/CMakeLists.txt
index 4a90cc4f475bfde3552794f60c25302e219ab800..cf383afbce8962b06b468a455652010a3f326da5 100644 (file)
--- a/source/lac/CMakeLists.txt
+++ b/source/lac/CMakeLists.txt
@@ -45,6 +45,7 @@ SET(_unity_include_src
   sparse_vanka.cc
   sparsity_pattern.cc
   sparsity_tools.cc
+  tensor_product_matrix.cc
   vector.cc
   vector_memory.cc
   )
@@ -86,6 +87,7 @@ SET(_inst
   solver.inst.in
   sparse_matrix_ez.inst.in
   sparse_matrix.inst.in
+  tensor_product_matrix.inst.in
   vector.inst.in
   vector_memory.inst.in
   vector_memory_release.inst.in

diff --git a/source/lac/tensor_product_matrix.cc b/source/lac/tensor_product_matrix.cc
new file mode 100644 (file)
index 0000000..223796f
--- /dev/null
+++ b/source/lac/tensor_product_matrix.cc
@@ -0,0 +1,29 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+#include <deal.II/lac/tensor_product_matrix.templates.h>
+
+DEAL_II_NAMESPACE_OPEN
+
+namespace internal
+{
+  namespace TensorProductMatrixSymmetricSum
+  {
+#include "tensor_product_matrix.inst"
+
+  }
+} // namespace internal
+
+DEAL_II_NAMESPACE_CLOSE

diff --git a/source/lac/tensor_product_matrix.inst.in b/source/lac/tensor_product_matrix.inst.in
new file mode 100644 (file)
index 0000000..f2f8446
--- /dev/null
+++ b/source/lac/tensor_product_matrix.inst.in
@@ -0,0 +1,58 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2022 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE.md at
+// the top level directory of deal.II.
+//
+// ---------------------------------------------------------------------
+
+
+for (deal_II_dimension : DIMENSIONS;
+     deal_II_scalar_vectorized : REAL_SCALARS_VECTORIZED)
+  {
+    template void select_vmult<1>(
+      deal_II_scalar_vectorized * dst,
+      const deal_II_scalar_vectorized *         src,
+      AlignedVector<deal_II_scalar_vectorized> &tmp,
+      const std::array<Table<2, deal_II_scalar_vectorized>, deal_II_dimension>
+        &mass_matrix,
+      const std::array<Table<2, deal_II_scalar_vectorized>, deal_II_dimension>
+        &derivative_matrix);
+
+    template void select_apply_inverse<1>(
+      deal_II_scalar_vectorized * dst,
+      const deal_II_scalar_vectorized *         src,
+      AlignedVector<deal_II_scalar_vectorized> &tmp,
+      const std::array<Table<2, deal_II_scalar_vectorized>, deal_II_dimension>
+        &                                  eigenvectors,
+      const std::array<AlignedVector<deal_II_scalar_vectorized>,
+                       deal_II_dimension> &eigenvalues);
+  }
+
+for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS)
+  {
+    template void select_vmult<1>(
+      deal_II_scalar * dst,
+      const deal_II_scalar *         src,
+      AlignedVector<deal_II_scalar> &tmp,
+      const std::array<Table<2, deal_II_scalar>, deal_II_dimension>
+        &mass_matrix,
+      const std::array<Table<2, deal_II_scalar>, deal_II_dimension>
+        &derivative_matrix);
+
+    template void select_apply_inverse<1>(
+      deal_II_scalar * dst,
+      const deal_II_scalar *         src,
+      AlignedVector<deal_II_scalar> &tmp,
+      const std::array<Table<2, deal_II_scalar>, deal_II_dimension>
+        &eigenvectors,
+      const std::array<AlignedVector<deal_II_scalar>, deal_II_dimension>
+        &eigenvalues);
+  }