From 8df7e16650757db75e1d76283dada003265739c2 Mon Sep 17 00:00:00 2001 From: Martin Kronbichler Date: Tue, 12 Oct 2021 15:09:40 +0200 Subject: [PATCH] Enable PreconditionChebyshev to update a sub-range of vector entries --- include/deal.II/lac/precondition.h | 230 ++++++++++++++++++++++++++--- 1 file changed, 207 insertions(+), 23 deletions(-) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index fdfa87d49c..2258633108 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -976,8 +976,8 @@ private: * *

Requirements on the templated classes

* - * The class MatrixType must be derived from Subscriptor because a - * SmartPointer to MatrixType is held in the class. In particular, this means + * The class `MatrixType` must be derived from Subscriptor because a + * SmartPointer to `MatrixType` is held in the class. In particular, this means * that the matrix object needs to persist during the lifetime of * PreconditionChebyshev. The preconditioner is held in a shared_ptr that is * copied into the AdditionalData member variable of the class, so the @@ -994,6 +994,47 @@ private: * entries that would be needed from the matrix alone), there is a backward * compatibility function that can extract the diagonal in case of a serial * computation. + * + *

Optimized operations with specific `MatrixType` argument

+ * + * This class enables to embed the vector updates into the matrix-vector + * product in case the `MatrixType` supports this. To this end, the + * `VectorType` needs to be of type LinearAlgebra::distributed::Vector, the + * `PreconditionerType` needs to be DiagonalMatrix, and the class `MatrixType` + * needs to provide a function with the signature + * @code + * void MatrixType::vmult( + * VectorType &, + * const VectorType &, + * const std::function &, + * const std::function &) const + * @endcode + * where the two given functions run before and after the matrix-vector + * product, respectively. They take as arguments a sub-range among the locally + * owned elements of the vector, and allow the matrix-vector product to return + * an index range that fulfills all the requirements. For the example of a + * class similar to the one in the step-37 tutorial program, the implementation + * is + * @code + * void + * vmult(LinearAlgebra::distributed::Vector & dst, + * const LinearAlgebra::distributed::Vector &src, + * const std::function + * &operation_before_loop, + * const std::function + * &operation_after_loop) const + * { + * data.cell_loop(&LaplaceOperator::local_apply, + * this, + * dst, + * src, + * operation_before_loop, + * operation_after_loop); + * } + * @endcode + * In terms of the Chebyshev iteration, the operation before the loop will + * set `dst` to zero, whereas the operation after the loop performs the + * iteration leading to $x^{n+1}$ described above. */ template , typename VectorType = Vector, @@ -2036,6 +2077,148 @@ namespace internal solution.swap(solution_old); } + // Detector class to find out whether we have a vmult function that takes + // two additional std::function objects, which we use to run the operation + // on slices of the vector during the matrix-vector product + template + struct supports_vmult_with_std_functions + { + private: + // this will work always + static bool + detect(...); + + // this detector will work only if we have + // "... MatrixType::vmult(VectorType, const VectorType, + // const std::function<...>&, const std::function<...>&) const" + template + static decltype(std::declval().vmult( + std::declval(), + std::declval(), + std::declval &>(), + std::declval &>())) + detect(const MatrixType2 &); + + public: + // finally here we check if our detector has void return type and + // fulfills additional requirements on the vector type and + // preconditioner. This will happen if the compiler can use the second + // detector, otherwise SFINAE let's it work with the more general first + // one that is bool + static const bool value = + !std::is_same()))>::value && + std::is_same>::value && + std::is_same< + VectorType, + LinearAlgebra::distributed::Vector>::value; + }; + + // We need to have a separate declaration for static const members + template + const bool supports_vmult_with_std_functions::value; + + template ::value, + int>::type * = nullptr> + inline void + vmult_and_update(const MatrixType & matrix, + const PreconditionerType &preconditioner, + const VectorType & rhs, + const unsigned int iteration_index, + const double factor1, + const double factor2, + VectorType & solution, + VectorType & solution_old, + VectorType & temp_vector1, + VectorType & temp_vector2) + { + if (iteration_index > 0) + matrix.vmult(temp_vector1, solution); + vector_updates(rhs, + preconditioner, + iteration_index, + factor1, + factor2, + solution_old, + temp_vector1, + temp_vector2, + solution); + } + + template ::value, + int>::type * = nullptr> + inline void + vmult_and_update(const MatrixType & matrix, + const PreconditionerType &preconditioner, + const VectorType & rhs, + const unsigned int iteration_index, + const double factor1, + const double factor2, + VectorType & solution, + VectorType & solution_old, + VectorType & temp_vector1, + VectorType &) + { + using Number = typename VectorType::value_type; + VectorUpdater updater(rhs.begin(), + preconditioner.get_vector().begin(), + iteration_index, + factor1, + factor2, + solution_old.begin(), + temp_vector1.begin(), + solution.begin()); + if (iteration_index > 0) + matrix.vmult( + temp_vector1, + solution, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'temp_vector1' before running the vmult + // operation + if (end_range > start_range) + std::memset(temp_vector1.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int start_range, const unsigned int end_range) { + if (end_range > start_range) + updater.apply_to_subrange(start_range, end_range); + }); + else + updater.apply_to_subrange(0U, solution.locally_owned_size()); + + // swap vectors x^{n+1}->x^{n}, given the updates in the function above + if (iteration_index == 0) + { + // nothing to do here because we can immediately write into the + // solution vector without remembering any of the other vectors + } + else if (iteration_index == 1) + { + solution.swap(temp_vector1); + solution_old.swap(temp_vector1); + } + else + solution.swap(solution_old); + } + template inline void initialize_preconditioner( @@ -2421,16 +2604,17 @@ PreconditionChebyshev::vmult( if (eigenvalues_are_initialized == false) estimate_eigenvalues(rhs); - internal::PreconditionChebyshevImplementation::vector_updates( - rhs, + internal::PreconditionChebyshevImplementation::vmult_and_update( + *matrix_ptr, *data.preconditioner, + rhs, 0, 0., 1. / theta, + solution, solution_old, temp_vector1, - temp_vector2, - solution); + temp_vector2); // if delta is zero, we do not need to iterate because the updates will be // zero @@ -2440,20 +2624,20 @@ PreconditionChebyshev::vmult( double rhok = delta / theta, sigma = theta / delta; for (unsigned int k = 0; k < data.degree - 1; ++k) { - matrix_ptr->vmult(temp_vector1, solution); const double rhokp = 1. / (2. * sigma - rhok); const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta; rhok = rhokp; - internal::PreconditionChebyshevImplementation::vector_updates( - rhs, + internal::PreconditionChebyshevImplementation::vmult_and_update( + *matrix_ptr, *data.preconditioner, + rhs, k + 1, factor1, factor2, + solution, solution_old, temp_vector1, - temp_vector2, - solution); + temp_vector2); } } @@ -2486,10 +2670,10 @@ PreconditionChebyshev::Tvmult( double rhok = delta / theta, sigma = theta / delta; for (unsigned int k = 0; k < data.degree - 1; ++k) { - matrix_ptr->Tvmult(temp_vector1, solution); const double rhokp = 1. / (2. * sigma - rhok); const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta; rhok = rhokp; + matrix_ptr->Tvmult(temp_vector1, solution); internal::PreconditionChebyshevImplementation::vector_updates( rhs, *data.preconditioner, @@ -2515,17 +2699,17 @@ PreconditionChebyshev::step( if (eigenvalues_are_initialized == false) estimate_eigenvalues(rhs); - matrix_ptr->vmult(temp_vector1, solution); - internal::PreconditionChebyshevImplementation::vector_updates( - rhs, + internal::PreconditionChebyshevImplementation::vmult_and_update( + *matrix_ptr, *data.preconditioner, + rhs, 1, 0., 1. / theta, + solution, solution_old, temp_vector1, - temp_vector2, - solution); + temp_vector2); if (data.degree < 2 || std::abs(delta) < 1e-40) return; @@ -2533,20 +2717,20 @@ PreconditionChebyshev::step( double rhok = delta / theta, sigma = theta / delta; for (unsigned int k = 0; k < data.degree - 1; ++k) { - matrix_ptr->vmult(temp_vector1, solution); const double rhokp = 1. / (2. * sigma - rhok); const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta; rhok = rhokp; - internal::PreconditionChebyshevImplementation::vector_updates( - rhs, + internal::PreconditionChebyshevImplementation::vmult_and_update( + *matrix_ptr, *data.preconditioner, + rhs, k + 2, factor1, factor2, + solution, solution_old, temp_vector1, - temp_vector2, - solution); + temp_vector2); } } @@ -2580,10 +2764,10 @@ PreconditionChebyshev::Tstep( double rhok = delta / theta, sigma = theta / delta; for (unsigned int k = 0; k < data.degree - 1; ++k) { - matrix_ptr->Tvmult(temp_vector1, solution); const double rhokp = 1. / (2. * sigma - rhok); const double factor1 = rhokp * rhok, factor2 = 2. * rhokp / delta; rhok = rhokp; + matrix_ptr->Tvmult(temp_vector1, solution); internal::PreconditionChebyshevImplementation::vector_updates( rhs, *data.preconditioner, -- 2.39.5