From 06186bdd6270f7bfbc2b4805f99b656778c0d053 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 31 May 2023 19:39:21 +0200 Subject: [PATCH] Don't use tmp vector --- include/deal.II/lac/tensor_product_matrix.h | 103 +++++------------- .../lac/tensor_product_matrix.templates.h | 27 +---- source/lac/tensor_product_matrix.inst.in | 6 +- tests/lac/tensor_product_matrix_08.cc | 13 +-- 4 files changed, 38 insertions(+), 111 deletions(-) diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 20550c3f00..57745af67d 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -209,28 +209,11 @@ public: * described in the main documentation of TensorProductMatrixSymmetricSum. * This function is operating on ArrayView to allow checks of * array bounds with respect to @p dst and @p src. - * - * @warning This function works on an internal temporal array, leading to - * increased memory consumption if many instances of this class are created, - * e.g., a different object on every cell with different underlying - * coefficients each. Furthermore, only one thread runs this function at a - * time (ensured internally with a mutex). If these two limitations are an - * issue, please consider the other version of this function. */ void apply_inverse(const ArrayView & dst, const ArrayView &src) const; - /** - * Same as above but letting the user provide a user-owned temporary array, - * resolving the two issues described above. This array is resized - * internally to the needed size. - */ - void - apply_inverse(const ArrayView & dst, - const ArrayView &src, - AlignedVector & tmp) const; - /** * Return the memory consumption of the allocated memory in this class. */ @@ -418,8 +401,7 @@ public: void apply_inverse(const unsigned int index, const ArrayView & dst_in, - const ArrayView &src_in, - AlignedVector & tmp_array) const; + const ArrayView &src_in) const; /** * Return the memory consumption of this class in bytes. @@ -819,10 +801,9 @@ namespace internal template void - apply_inverse(Number * dst, - const Number * src, - AlignedVector &tmp, - const unsigned int n_rows_1d_non_templated, + apply_inverse(Number * dst, + const Number * src, + const unsigned int n_rows_1d_non_templated, const std::array &eigenvectors, const std::array &eigenvalues, const Number *inverted_eigenvalues = nullptr) @@ -830,10 +811,6 @@ namespace internal const unsigned int n_rows_1d = n_rows_1d_templated == 0 ? n_rows_1d_non_templated : n_rows_1d_templated; - const unsigned int n = Utilities::fixed_power(n_rows_1d); - - tmp.resize_fast(n); - Number *t = tmp.begin(); internal::EvaluatorTensorProduct(S, src, t); + eval.template apply<0, true, false>(S, src, dst); for (unsigned int i = 0; i < n_rows_1d; ++i) if (inverted_eigenvalues) - t[i] *= inverted_eigenvalues[i]; + dst[i] *= inverted_eigenvalues[i]; else - t[i] /= eigenvalues[0][i]; + dst[i] /= eigenvalues[0][i]; - eval.template apply<0, false, false>(S, t, dst); + eval.template apply<0, false, false>(S, dst, dst); } else if (dim == 2) { const Number *S0 = eigenvectors[0]; const Number *S1 = eigenvectors[1]; - eval.template apply<0, true, false>(S0, src, t); - eval.template apply<1, true, false>(S1, t, dst); + eval.template apply<0, true, false>(S0, src, dst); + eval.template apply<1, true, false>(S1, dst, dst); for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1) for (unsigned int i0 = 0; i0 < n_rows_1d; ++i0, ++c) @@ -875,8 +852,8 @@ namespace internal else 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); + eval.template apply<1, false, false>(S1, dst, dst); + eval.template apply<0, false, false>(S0, dst, dst); } else if (dim == 3) @@ -884,22 +861,22 @@ namespace internal const Number *S0 = eigenvectors[0]; const Number *S1 = eigenvectors[1]; const Number *S2 = eigenvectors[2]; - 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); + eval.template apply<0, true, false>(S0, src, dst); + eval.template apply<1, true, false>(S1, dst, dst); + eval.template apply<2, true, false>(S2, dst, dst); 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) if (inverted_eigenvalues) - t[c] *= inverted_eigenvalues[c]; + dst[c] *= inverted_eigenvalues[c]; else - t[c] /= (eigenvalues[2][i2] + eigenvalues[1][i1] + - eigenvalues[0][i0]); + dst[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); + eval.template apply<2, false, false>(S2, dst, dst); + eval.template apply<1, false, false>(S1, dst, dst); + eval.template apply<0, false, false>(S0, dst, dst); } else @@ -923,7 +900,6 @@ namespace internal void select_apply_inverse(Number * dst, const Number * src, - AlignedVector & tmp, const unsigned int n_rows_1d, const std::array &eigenvectors, const std::array &eigenvalues, @@ -1016,19 +992,6 @@ inline void TensorProductMatrixSymmetricSum::apply_inverse( const ArrayView & dst_view, const ArrayView &src_view) const -{ - std::lock_guard lock(this->mutex); - this->apply_inverse(dst_view, src_view, this->tmp_array); -} - - - -template -inline void -TensorProductMatrixSymmetricSum::apply_inverse( - const ArrayView & dst_view, - const ArrayView &src_view, - AlignedVector & tmp_array) const { AssertDimension(dst_view.size(), this->n()); AssertDimension(src_view.size(), this->m()); @@ -1049,10 +1012,10 @@ TensorProductMatrixSymmetricSum::apply_inverse( if (n_rows_1d != -1) internal::TensorProductMatrixSymmetricSum::apply_inverse< n_rows_1d == -1 ? 0 : n_rows_1d>( - dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); + dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues); else internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( - dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); + dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues); } @@ -1518,8 +1481,7 @@ void TensorProductMatrixSymmetricSumCollection:: apply_inverse(const unsigned int index, const ArrayView & dst_in, - const ArrayView &src_in, - AlignedVector & tmp_array) const + const ArrayView &src_in) const { Number * dst = dst_in.begin(); const Number *src = src_in.begin(); @@ -1545,20 +1507,11 @@ TensorProductMatrixSymmetricSumCollection:: if (n_rows_1d != -1) internal::TensorProductMatrixSymmetricSum::apply_inverse< - n_rows_1d == -1 ? 0 : n_rows_1d>(dst, - src, - tmp_array, - n_rows_1d_non_templated, - eigenvectors, - eigenvalues); + n_rows_1d == -1 ? 0 : n_rows_1d>( + dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues); else internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( - dst, - src, - tmp_array, - n_rows_1d_non_templated, - eigenvectors, - eigenvalues); + dst, src, n_rows_1d_non_templated, eigenvectors, eigenvalues); } else { @@ -1595,7 +1548,6 @@ TensorProductMatrixSymmetricSumCollection:: internal::TensorProductMatrixSymmetricSum::apply_inverse< n_rows_1d == -1 ? 0 : n_rows_1d>(dst, src, - tmp_array, n_rows_1d_non_templated, eigenvectors, {}, @@ -1604,7 +1556,6 @@ TensorProductMatrixSymmetricSumCollection:: internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( dst, src, - tmp_array, n_rows_1d_non_templated, eigenvectors, {}, diff --git a/include/deal.II/lac/tensor_product_matrix.templates.h b/include/deal.II/lac/tensor_product_matrix.templates.h index 90a16f2e32..21a52a6e84 100644 --- a/include/deal.II/lac/tensor_product_matrix.templates.h +++ b/include/deal.II/lac/tensor_product_matrix.templates.h @@ -66,37 +66,20 @@ namespace internal void select_apply_inverse(Number * dst, const Number * src, - AlignedVector & tmp, const unsigned int n_rows_1d, const std::array &eigenvectors, const std::array &eigenvalues, const Number *inverted_eigenvalues) { if (n_rows_1d_templated == n_rows_1d) - apply_inverse(dst, - src, - tmp, - n_rows_1d, - eigenvectors, - eigenvalues, - inverted_eigenvalues); + apply_inverse( + dst, src, n_rows_1d, eigenvectors, eigenvalues, inverted_eigenvalues); else if (n_rows_1d_templated < FDM_N_ROWS_MAX) select_apply_inverse( - dst, - src, - tmp, - n_rows_1d, - eigenvectors, - eigenvalues, - inverted_eigenvalues); + dst, src, n_rows_1d, eigenvectors, eigenvalues, inverted_eigenvalues); else - apply_inverse<0>(dst, - src, - tmp, - n_rows_1d, - eigenvectors, - eigenvalues, - inverted_eigenvalues); + apply_inverse<0>( + dst, src, n_rows_1d, eigenvectors, eigenvalues, inverted_eigenvalues); } } // namespace TensorProductMatrixSymmetricSum } // namespace internal diff --git a/source/lac/tensor_product_matrix.inst.in b/source/lac/tensor_product_matrix.inst.in index ec5c1c41ca..a7e8b3dd3d 100644 --- a/source/lac/tensor_product_matrix.inst.in +++ b/source/lac/tensor_product_matrix.inst.in @@ -29,9 +29,8 @@ for (deal_II_dimension : DIMENSIONS; template void select_apply_inverse<1>( deal_II_scalar_vectorized * dst, - const deal_II_scalar_vectorized * src, - AlignedVector &tmp, - const unsigned int n_rows, + const deal_II_scalar_vectorized *src, + const unsigned int n_rows, const std::array &eigenvectors, const std::array @@ -53,7 +52,6 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar : REAL_SCALARS) template void select_apply_inverse<1>( deal_II_scalar * dst, const deal_II_scalar * src, - AlignedVector & tmp, const unsigned int n_rows, const std::array &eigenvectors, const std::array &eigenvalues, diff --git a/tests/lac/tensor_product_matrix_08.cc b/tests/lac/tensor_product_matrix_08.cc index 9ee0b7b415..481fc7afc9 100644 --- a/tests/lac/tensor_product_matrix_08.cc +++ b/tests/lac/tensor_product_matrix_08.cc @@ -130,7 +130,6 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) AlignedVector src(fe.n_dofs_per_cell()); AlignedVector dst(fe.n_dofs_per_cell()); - AlignedVector tmp; Table<2, Number> matrix_0(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); Table<2, Number> matrix_1(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); Table<2, Number> matrix_2(fe.n_dofs_per_cell(), fe.n_dofs_per_cell()); @@ -145,29 +144,25 @@ do_test_mesh(const Mapping &mapping, const Triangulation &tria) collection_0.apply_inverse(cell, make_array_view(dst.begin(), dst.end()), - make_array_view(src.begin(), src.end()), - tmp); + make_array_view(src.begin(), src.end())); for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) matrix_0[j][i] = dst[j]; collection_1.apply_inverse(cell, make_array_view(dst.begin(), dst.end()), - make_array_view(src.begin(), src.end()), - tmp); + make_array_view(src.begin(), src.end())); for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) matrix_1[j][i] = dst[j]; collection_2.apply_inverse(cell, make_array_view(dst.begin(), dst.end()), - make_array_view(src.begin(), src.end()), - tmp); + make_array_view(src.begin(), src.end())); for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) matrix_2[j][i] = dst[j]; collection_3.apply_inverse(cell, make_array_view(dst.begin(), dst.end()), - make_array_view(src.begin(), src.end()), - tmp); + make_array_view(src.begin(), src.end())); for (unsigned int j = 0; j < fe.n_dofs_per_cell(); ++j) matrix_3[j][i] = dst[j]; } -- 2.39.5