From 162c0d038bc728916c150d0fce73d073118cd3ed Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 28 Aug 2022 17:52:59 +0200 Subject: [PATCH] Precompile TensorProductMatrixSymmetricSum kernels --- include/deal.II/lac/tensor_product_matrix.h | 291 +++++++++++------- .../lac/tensor_product_matrix.templates.h | 86 ++++++ source/lac/CMakeLists.txt | 2 + source/lac/tensor_product_matrix.cc | 29 ++ source/lac/tensor_product_matrix.inst.in | 58 ++++ 5 files changed, 353 insertions(+), 113 deletions(-) create mode 100644 include/deal.II/lac/tensor_product_matrix.templates.h create mode 100644 source/lac/tensor_product_matrix.cc create mode 100644 source/lac/tensor_product_matrix.inst.in diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index ffaa3de04c..0d2ef2a199 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -21,6 +21,7 @@ #include #include +#include #include @@ -502,6 +503,167 @@ TensorProductMatrixSymmetricSumBase::n() const } +namespace internal +{ + namespace TensorProductMatrixSymmetricSum + { + template + void + vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) + { + const unsigned int n_rows_1d = mass_matrix[0].n_rows(); + const unsigned int n = Utilities::fixed_power(n_rows_1d); + + tmp.resize_fast(n * 2); + Number *t = tmp.begin(); + + internal::EvaluatorTensorProduct + 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 + void + apply_inverse(Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> & eigenvectors, + const std::array, dim> &eigenvalues) + { + const unsigned int n_rows_1d = eigenvectors[0].n_rows(); + const unsigned int n = Utilities::fixed_power(n_rows_1d); + + tmp.resize_fast(n); + Number *t = tmp.begin(); + + internal::EvaluatorTensorProduct + 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 + void + select_vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix); + + + + template + void + select_apply_inverse( + Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> & eigenvectors, + const std::array, dim> &eigenvalues); + } // namespace TensorProductMatrixSymmetricSum + +} // namespace internal + + template inline void @@ -512,61 +674,17 @@ TensorProductMatrixSymmetricSumBase::vmult( AssertDimension(dst_view.size(), this->m()); AssertDimension(src_view.size(), this->n()); std::lock_guard lock(this->mutex); - const unsigned int n = Utilities::fixed_power( - 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 - eval(AlignedVector{}, - AlignedVector{}, - AlignedVector{}, - 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::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(n)); - constexpr int kernel_size = n_rows_1d > 0 ? n_rows_1d : 0; - internal::EvaluatorTensorProduct - eval(AlignedVector(), - AlignedVector(), - AlignedVector(), - 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 index 0000000000..4bbb26fec0 --- /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 + +#include + +#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 + void + select_vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> &mass_matrix, + const std::array, dim> &derivative_matrix) + { + const int n_rows_1d = mass_matrix[0].n_rows(); + + if (n_rows_1d_templated == n_rows_1d) + vmult( + dst, src, tmp, mass_matrix, derivative_matrix); + else if (n_rows_1d_templated < FDM_DEGREE_MAX) + select_vmult( + dst, src, tmp, mass_matrix, derivative_matrix); + else + vmult<0>(dst, src, tmp, mass_matrix, derivative_matrix); + } + + + + template + void + select_apply_inverse( + Number * dst, + const Number * src, + AlignedVector & tmp, + const std::array, dim> & eigenvectors, + const std::array, dim> &eigenvalues) + { + const int n_rows_1d = eigenvectors[0].n_rows(); + + if (n_rows_1d_templated == n_rows_1d) + apply_inverse( + dst, src, tmp, eigenvectors, eigenvalues); + else if (n_rows_1d_templated < FDM_DEGREE_MAX) + select_apply_inverse( + 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 4a90cc4f47..cf383afbce 100644 --- 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 index 0000000000..223796f8a3 --- /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_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 index 0000000000..f2f84467d8 --- /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 &tmp, + const std::array, deal_II_dimension> + &mass_matrix, + const std::array, deal_II_dimension> + &derivative_matrix); + + template void select_apply_inverse<1>( + deal_II_scalar_vectorized * dst, + const deal_II_scalar_vectorized * src, + AlignedVector &tmp, + const std::array, deal_II_dimension> + & eigenvectors, + const std::array, + 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 &tmp, + const std::array, deal_II_dimension> + &mass_matrix, + const std::array, deal_II_dimension> + &derivative_matrix); + + template void select_apply_inverse<1>( + deal_II_scalar * dst, + const deal_II_scalar * src, + AlignedVector &tmp, + const std::array, deal_II_dimension> + &eigenvectors, + const std::array, deal_II_dimension> + &eigenvalues); + } -- 2.39.5