From 76fc8de6f29385df6b60b97babf8ff0c86271022 Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Wed, 7 Sep 2022 12:04:49 +0200 Subject: [PATCH] DiagonalMatrix: specialize for LA:d:V --- include/deal.II/lac/diagonal_matrix.h | 52 ++++++++++++++++- include/deal.II/lac/precondition.h | 83 ++++++++++++++------------- 2 files changed, 93 insertions(+), 42 deletions(-) diff --git a/include/deal.II/lac/diagonal_matrix.h b/include/deal.II/lac/diagonal_matrix.h index 0fd53c312f..4c8a0c75ac 100644 --- a/include/deal.II/lac/diagonal_matrix.h +++ b/include/deal.II/lac/diagonal_matrix.h @@ -24,6 +24,20 @@ DEAL_II_NAMESPACE_OPEN +// forward declarations +#ifndef DOXYGEN +template +class Vector; +namespace LinearAlgebra +{ + namespace distributed + { + template + class Vector; + } // namespace distributed +} // namespace LinearAlgebra +#endif + /** * This class represents a n x n diagonal matrix based on a vector of * size n. The matrix-vector products are realized by @p @@ -386,12 +400,46 @@ DiagonalMatrix::add(const size_type i, +namespace internal +{ + namespace DiagonalMatrix + { + template + void + assign_and_scale(VectorType & dst, + const VectorType &src, + const VectorType &diagonal) + { + dst = src; + dst.scale(diagonal); + } + + template + void + assign_and_scale( + LinearAlgebra::distributed::Vector & dst, + const LinearAlgebra::distributed::Vector &src, + const LinearAlgebra::distributed::Vector + &diagonal) + { + const auto dst_ptr = dst.begin(); + const auto src_ptr = src.begin(); + const auto diagonal_ptr = diagonal.begin(); + + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < src.locally_owned_size(); ++i) + dst_ptr[i] = src_ptr[i] * diagonal_ptr[i]; + } + } // namespace DiagonalMatrix +} // namespace internal + + + template void DiagonalMatrix::vmult(VectorType &dst, const VectorType &src) const { - dst = src; - dst.scale(diagonal); + internal::DiagonalMatrix::assign_and_scale(dst, src, diagonal); } diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 2972ddc408..b5702d9959 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -542,7 +542,8 @@ namespace internal typename PreconditionerType> constexpr bool has_vmult_with_std_functions = is_supported_operation && - std::is_same>::value && + std::is_same>::value && (std::is_same>::value || std::is_same< @@ -1175,19 +1176,19 @@ namespace internal // 2) specialized implementation for inverse-diagonal preconditioner template ::value && - !has_vmult_with_std_functions>, - VectorType> * = nullptr> + std::enable_if_t::value && + !has_vmult_with_std_functions< + MatrixType, + VectorType, + dealii::DiagonalMatrix>, + VectorType> * = nullptr> void - step_operations(const MatrixType & A, - const DiagonalMatrix &preconditioner, - VectorType & dst, - const VectorType & src, - const double relaxation, - VectorType & tmp, + step_operations(const MatrixType & A, + const dealii::DiagonalMatrix &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & tmp, VectorType &, const unsigned int i, const bool transposed) @@ -1247,19 +1248,19 @@ namespace internal // matrix that accepts ranges template ::value && - has_vmult_with_std_functions>, - VectorType> * = nullptr> + std::enable_if_t::value && + has_vmult_with_std_functions< + MatrixType, + VectorType, + dealii::DiagonalMatrix>, + VectorType> * = nullptr> void - step_operations(const MatrixType & A, - const DiagonalMatrix &preconditioner, - VectorType & dst, - const VectorType & src, - const double relaxation, - VectorType & tmp, + step_operations(const MatrixType & A, + const dealii::DiagonalMatrix &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & tmp, VectorType &, const unsigned int i, const bool transposed) @@ -2871,15 +2872,16 @@ namespace internal // selection for diagonal matrix around deal.II vector template inline void - vector_updates(const ::dealii::Vector & rhs, - const DiagonalMatrix<::dealii::Vector> &jacobi, - const unsigned int iteration_index, - const double factor1, - const double factor2, - ::dealii::Vector &solution_old, - ::dealii::Vector &temp_vector1, - ::dealii::Vector &, - ::dealii::Vector &solution) + vector_updates( + const ::dealii::Vector & rhs, + const dealii::DiagonalMatrix<::dealii::Vector> &jacobi, + const unsigned int iteration_index, + const double factor1, + const double factor2, + ::dealii::Vector & solution_old, + ::dealii::Vector & temp_vector1, + ::dealii::Vector &, + ::dealii::Vector &solution) { VectorUpdater upd(rhs.begin(), jacobi.get_vector().begin(), @@ -2909,7 +2911,7 @@ namespace internal inline void vector_updates( const LinearAlgebra::distributed::Vector &rhs, - const DiagonalMatrix< + const dealii::DiagonalMatrix< LinearAlgebra::distributed::Vector> &jacobi, const unsigned int iteration_index, const double factor1, @@ -3053,13 +3055,14 @@ namespace internal template inline void initialize_preconditioner( - const MatrixType & matrix, - std::shared_ptr> &preconditioner) + const MatrixType & matrix, + std::shared_ptr> &preconditioner) { if (preconditioner.get() == nullptr || preconditioner->m() != matrix.m()) { if (preconditioner.get() == nullptr) - preconditioner = std::make_shared>(); + preconditioner = + std::make_shared>(); Assert( preconditioner->m() == 0, @@ -3438,8 +3441,8 @@ PreconditionChebyshev:: // We do not need the second temporary vector in case we have a // DiagonalMatrix as preconditioner and use deal.II's own vectors using NumberType = typename VectorType::value_type; - if (std::is_same>::value == - false || + if (std::is_same>::value == false || (std::is_same>::value == false && ((std::is_same