From: Peter Munch Date: Mon, 19 Sep 2022 11:21:12 +0000 (+0200) Subject: TensorProductMatrixSymmetricSum: switch to pointers X-Git-Tag: v9.5.0-rc1~952^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14288%2Fhead;p=dealii.git TensorProductMatrixSymmetricSum: switch to pointers --- diff --git a/include/deal.II/lac/tensor_product_matrix.h b/include/deal.II/lac/tensor_product_matrix.h index 8734f694a2..850e3249fc 100644 --- a/include/deal.II/lac/tensor_product_matrix.h +++ b/include/deal.II/lac/tensor_product_matrix.h @@ -182,10 +182,27 @@ 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 vmult(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 + vmult(const ArrayView & dst, + const ArrayView &src, + AlignedVector & tmp) const; + /** * Implements a matrix-vector product with the underlying matrix as * described in the main documentation of TensorProductMatrixSymmetricSum. @@ -195,16 +212,16 @@ public: * @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 run this function at once - * (ensured internally with a mutex). If these two limitations are an issue, - * please consider the other version of this function. + * 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 the user can provide a user-owned temporal array, + * 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. */ @@ -489,13 +506,16 @@ namespace internal template void - vmult(Number * dst, - const Number * src, - AlignedVector & tmp, - const std::array, dim> &mass_matrix, - const std::array, dim> &derivative_matrix) + vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const unsigned int n_rows_1d_non_templated, + const std::array &mass_matrix, + const std::array &derivative_matrix) { - const unsigned int n_rows_1d = mass_matrix[0].n_rows(); + 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 * 2); @@ -510,16 +530,16 @@ namespace internal if (dim == 1) { - const Number *A = &derivative_matrix[0](0, 0); + const Number *A = derivative_matrix[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); + const Number *A0 = derivative_matrix[0]; + const Number *M0 = mass_matrix[0]; + const Number *A1 = derivative_matrix[1]; + const Number *M1 = mass_matrix[1]; 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); @@ -528,12 +548,12 @@ namespace internal 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); + const Number *A0 = derivative_matrix[0]; + const Number *M0 = mass_matrix[0]; + const Number *A1 = derivative_matrix[1]; + const Number *M1 = mass_matrix[1]; + const Number *A2 = derivative_matrix[2]; + const Number *M2 = mass_matrix[2]; 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); @@ -551,13 +571,16 @@ namespace internal template void - apply_inverse(Number * dst, - const Number * src, - AlignedVector & tmp, - const std::array, dim> & eigenvectors, - const std::array, dim> &eigenvalues) + apply_inverse(Number * dst, + const Number * src, + AlignedVector &tmp, + const unsigned int n_rows_1d_non_templated, + const std::array &eigenvectors, + const std::array &eigenvalues) { - const unsigned int n_rows_1d = eigenvectors[0].n_rows(); + 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); @@ -577,7 +600,7 @@ namespace internal // rows correspond to dofs whereas columns to eigenvalue indices! if (dim == 1) { - const Number *S = &eigenvectors[0](0, 0); + const Number *S = eigenvectors[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]; @@ -586,8 +609,8 @@ namespace internal else if (dim == 2) { - const Number *S0 = &(eigenvectors[0](0, 0)); - const Number *S1 = &(eigenvectors[1](0, 0)); + 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); for (unsigned int i1 = 0, c = 0; i1 < n_rows_1d; ++i1) @@ -599,9 +622,9 @@ namespace internal 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); + 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); @@ -623,22 +646,23 @@ namespace internal template void - select_vmult(Number * dst, - const Number * src, - AlignedVector & tmp, - const std::array, dim> &mass_matrix, - const std::array, dim> &derivative_matrix); + select_vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const unsigned int n_rows_1d, + const std::array &mass_matrix, + const std::array &derivative_matrix); template void - select_apply_inverse( - Number * dst, - const Number * src, - AlignedVector & tmp, - const std::array, dim> & eigenvectors, - const std::array, dim> &eigenvalues); + select_apply_inverse(Number * dst, + const Number * src, + AlignedVector & tmp, + const unsigned int n_rows_1d, + const std::array &eigenvectors, + const std::array &eigenvalues); } // namespace TensorProductMatrixSymmetricSum } // namespace internal @@ -672,21 +696,52 @@ inline void TensorProductMatrixSymmetricSum::vmult( const ArrayView & dst_view, const ArrayView &src_view) const +{ + std::lock_guard lock(this->mutex); + this->vmult(dst_view, src_view, this->tmp_array); +} + + + +template +inline void +TensorProductMatrixSymmetricSum::vmult( + const ArrayView & dst_view, + const ArrayView &src_view, + AlignedVector & tmp_array) const { AssertDimension(dst_view.size(), this->m()); AssertDimension(src_view.size(), this->n()); - std::lock_guard lock(this->mutex); Number * dst = dst_view.begin(); const Number *src = src_view.begin(); + std::array mass_matrix, derivative_matrix; + + for (unsigned int d = 0; d < dim; ++d) + { + mass_matrix[d] = &this->mass_matrix[d](0, 0); + derivative_matrix[d] = &this->derivative_matrix[d](0, 0); + } + + const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows(); + if (n_rows_1d != -1) internal::TensorProductMatrixSymmetricSum::vmult< - n_rows_1d == -1 ? 0 : n_rows_1d>( - dst, src, tmp_array, mass_matrix, derivative_matrix); + n_rows_1d == -1 ? 0 : n_rows_1d>(dst, + src, + tmp_array, + n_rows_1d_non_templated, + mass_matrix, + derivative_matrix); else internal::TensorProductMatrixSymmetricSum::select_vmult<1>( - dst, src, tmp_array, mass_matrix, derivative_matrix); + dst, + src, + tmp_array, + n_rows_1d_non_templated, + mass_matrix, + derivative_matrix); } @@ -716,13 +771,23 @@ TensorProductMatrixSymmetricSum::apply_inverse( Number * dst = dst_view.begin(); const Number *src = src_view.begin(); + std::array eigenvectors, eigenvalues; + + for (unsigned int d = 0; d < dim; ++d) + { + eigenvectors[d] = &this->eigenvectors[d](0, 0); + eigenvalues[d] = this->eigenvalues[d].data(); + } + + const unsigned int n_rows_1d_non_templated = this->mass_matrix[0].n_rows(); + if (n_rows_1d != -1) internal::TensorProductMatrixSymmetricSum::apply_inverse< n_rows_1d == -1 ? 0 : n_rows_1d>( - dst, src, tmp_array, eigenvectors, eigenvalues); + dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); else internal::TensorProductMatrixSymmetricSum::select_apply_inverse<1>( - dst, src, tmp_array, eigenvectors, eigenvalues); + dst, src, tmp_array, n_rows_1d_non_templated, eigenvectors, eigenvalues); } diff --git a/include/deal.II/lac/tensor_product_matrix.templates.h b/include/deal.II/lac/tensor_product_matrix.templates.h index 4bbb26fec0..02cd64cfa7 100644 --- a/include/deal.II/lac/tensor_product_matrix.templates.h +++ b/include/deal.II/lac/tensor_product_matrix.templates.h @@ -37,45 +37,42 @@ namespace internal { template void - select_vmult(Number * dst, - const Number * src, - AlignedVector & tmp, - const std::array, dim> &mass_matrix, - const std::array, dim> &derivative_matrix) + select_vmult(Number * dst, + const Number * src, + AlignedVector & tmp, + const unsigned int n_rows_1d, + const std::array &mass_matrix, + const std::array &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); + dst, src, tmp, n_rows_1d, mass_matrix, derivative_matrix); else if (n_rows_1d_templated < FDM_DEGREE_MAX) select_vmult( - dst, src, tmp, mass_matrix, derivative_matrix); + dst, src, tmp, n_rows_1d, mass_matrix, derivative_matrix); else - vmult<0>(dst, src, tmp, mass_matrix, derivative_matrix); + vmult<0>(dst, src, tmp, n_rows_1d, 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) + 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 int n_rows_1d = eigenvectors[0].n_rows(); - if (n_rows_1d_templated == n_rows_1d) apply_inverse( - dst, src, tmp, eigenvectors, eigenvalues); + dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); else if (n_rows_1d_templated < FDM_DEGREE_MAX) select_apply_inverse( - dst, src, tmp, eigenvectors, eigenvalues); + dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); else - apply_inverse<0>(dst, src, tmp, eigenvectors, eigenvalues); + apply_inverse<0>(dst, src, tmp, n_rows_1d, eigenvectors, eigenvalues); } } // namespace TensorProductMatrixSymmetricSum } // namespace internal diff --git a/source/lac/tensor_product_matrix.inst.in b/source/lac/tensor_product_matrix.inst.in index f2f84467d8..6607de8aca 100644 --- a/source/lac/tensor_product_matrix.inst.in +++ b/source/lac/tensor_product_matrix.inst.in @@ -21,38 +21,39 @@ for (deal_II_dimension : DIMENSIONS; deal_II_scalar_vectorized * dst, const deal_II_scalar_vectorized * src, AlignedVector &tmp, - const std::array, deal_II_dimension> + const unsigned int n_rows, + const std::array &mass_matrix, - const std::array, deal_II_dimension> + const std::array &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); + const unsigned int n_rows, + const std::array + &eigenvectors, + const std::array + &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> + const deal_II_scalar * src, + AlignedVector & tmp, + const unsigned int n_rows, + const std::array &mass_matrix, + const std::array &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); + const deal_II_scalar * src, + AlignedVector & tmp, + const unsigned int n_rows, + const std::array &eigenvectors, + const std::array &eigenvalues); }