From: Peter Munch Date: Sun, 21 Aug 2022 20:24:52 +0000 (+0200) Subject: Optimize PreconditionRelaxation X-Git-Tag: v9.5.0-rc1~1019^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c7f1f49452a9cea40eadd80aa57a3a7a49f259b4;p=dealii.git Optimize PreconditionRelaxation --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index cc428ae8db..c75ede1e54 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -525,6 +525,31 @@ protected: namespace internal { + // a helper type-trait that leverage SFINAE to figure out if MatrixType has + // ... MatrixType::vmult(VectorType &, const VectorType&, + // std::function<...>, std::function<...>) const + template + using vmult_functions_t = decltype(std::declval().vmult( + std::declval(), + std::declval(), + std::declval< + const std::function &>(), + std::declval< + const std::function &>())); + + template + constexpr bool has_vmult_with_std_functions = + is_supported_operation && + std::is_same>::value && + (std::is_same>::value || + std::is_same< + VectorType, + LinearAlgebra::distributed::Vector>::value); + namespace PreconditionRelaxation { template @@ -1051,45 +1076,157 @@ namespace internal template ::value, VectorType> - * = nullptr> + std::enable_if_t< + !IsBlockVector::value && + !has_vmult_with_std_functions>, + VectorType> * = nullptr> void step_operations(const MatrixType & A, const DiagonalMatrix &preconditioner, VectorType & dst, const VectorType & src, const double relaxation, - VectorType & residual, + VectorType & tmp, VectorType &, const unsigned int i, const bool transposed) { + using Number = typename VectorType::value_type; + if (i == 0) { - const auto dst_ptr = dst.begin(); - const auto src_ptr = src.begin(); - const auto diag_ptr = preconditioner.get_vector().begin(); + Number * dst_ptr = dst.begin(); + const Number *src_ptr = src.begin(); + const Number *diag_ptr = preconditioner.get_vector().begin(); - for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) - dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i]; + if (relaxation == 1.0) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] = src_ptr[i] * diag_ptr[i]; + } + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i]; + } } else { - residual.reinit(src, true); + tmp.reinit(src, true); - const auto dst_ptr = dst.begin(); - const auto src_ptr = src.begin(); - const auto residual_ptr = residual.begin(); - const auto diag_ptr = preconditioner.get_vector().begin(); + Number * dst_ptr = dst.begin(); + const Number *src_ptr = src.begin(); + const Number *tmp_ptr = tmp.begin(); + const Number *diag_ptr = preconditioner.get_vector().begin(); if (transposed) - A.Tvmult(residual, dst); + Tvmult(A, tmp, dst); else - A.vmult(residual, dst); + A.vmult(tmp, dst); - for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) - dst_ptr[i] += - relaxation * (src_ptr[i] - residual_ptr[i]) * diag_ptr[i]; + if (relaxation == 1.0) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] += (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i]; + } + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] += + relaxation * (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i]; + } + } + } + + template ::value && + has_vmult_with_std_functions>, + VectorType> * = nullptr> + void + step_operations(const MatrixType & A, + const DiagonalMatrix &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & tmp, + VectorType &, + const unsigned int i, + const bool transposed) + { + using Number = typename VectorType::value_type; + + if (i == 0) + { + Number * dst_ptr = dst.begin(); + const Number *src_ptr = src.begin(); + const Number *diag_ptr = preconditioner.get_vector().begin(); + + if (relaxation == 1.0) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] = src_ptr[i] * diag_ptr[i]; + } + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (unsigned int i = 0; i < dst.locally_owned_size(); ++i) + dst_ptr[i] = relaxation * src_ptr[i] * diag_ptr[i]; + } + } + else + { + tmp.reinit(src, true); + + Assert(transposed == false, ExcNotImplemented()); + + A.vmult( + tmp, + dst, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'tmp' before running the vmult + // operation + if (end_range > start_range) + std::memset(tmp.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int begin, const unsigned int end) { + const Number *dst_ptr = dst.begin(); + const Number *src_ptr = src.begin(); + Number * tmp_ptr = tmp.begin(); + const Number *diag_ptr = preconditioner.get_vector().begin(); + + // for efficiency reason, write back to temp_vector that is + // already read (avoid read-for-ownership) + if (relaxation == 1.0) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = begin; i < end; ++i) + tmp_ptr[i] = + dst_ptr[i] + (src_ptr[i] - tmp_ptr[i]) * diag_ptr[i]; + } + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = begin; i < end; ++i) + tmp_ptr[i] = dst_ptr[i] + relaxation * + (src_ptr[i] - tmp_ptr[i]) * + diag_ptr[i]; + } + }); + + tmp.swap(dst); } } @@ -2536,29 +2673,6 @@ namespace internal solution.swap(solution_old); } - // a helper type-trait that leverage SFINAE to figure out if MatrixType has - // ... MatrixType::vmult(VectorType &, const VectorType&, - // std::function<...>, std::function<...>) const - template - using vmult_functions_t = decltype(std::declval().vmult( - std::declval(), - std::declval(), - std::declval< - const std::function &>(), - std::declval &>())); - - template - constexpr bool has_vmult_with_std_functions = - is_supported_operation - && std::is_same>::value - &&std::is_same< - VectorType, - LinearAlgebra::distributed::Vector>::value; - // We need to have a separate declaration for static const members template < diff --git a/tests/lac/precondition_relaxation_01.cc b/tests/lac/precondition_relaxation_01.cc index 2127481b66..884a4779a2 100644 --- a/tests/lac/precondition_relaxation_01.cc +++ b/tests/lac/precondition_relaxation_01.cc @@ -63,25 +63,64 @@ private: DiagonalMatrix diagonal_matrix; }; +template +class MySparseMatrix : public Subscriptor +{ +public: + MySparseMatrix(const SparseMatrixType &sparse_matrix) + : sparse_matrix(sparse_matrix) + {} + + template + void + vmult(VectorType & dst, + const VectorType &src, + const std::function + &operation_before_matrix_vector_product, + const std::function + &operation_after_matrix_vector_product) const + { + operation_before_matrix_vector_product(0, src.size()); + + sparse_matrix.vmult(dst, src); + + operation_after_matrix_vector_product(0, src.size()); + } + +private: + const SparseMatrixType &sparse_matrix; +}; template std::tuple -test(const PreconditionerType &preconditioner, const VectorType &src) +test(const PreconditionerType &preconditioner, + const VectorType & src, + const bool test_transposed = true) { VectorType dst; dst.reinit(src); + dst = 1.0; preconditioner.vmult(dst, src); const double norm_0 = dst.l2_norm(); + if (test_transposed) + { + dst = 1.0; + preconditioner.Tvmult(dst, src); + } + const double norm_2 = dst.l2_norm(); + + dst = 1.0; preconditioner.step(dst, src); const double norm_1 = dst.l2_norm(); - preconditioner.Tvmult(dst, src); - const double norm_2 = dst.l2_norm(); - - preconditioner.Tstep(dst, src); + if (test_transposed) + { + dst = 1.0; + preconditioner.Tstep(dst, src); + } const double norm_3 = dst.l2_norm(); return std::tuple{norm_0, @@ -159,6 +198,30 @@ main() results.emplace_back(test(preconditioner, src)); } + { + // Test PreconditionRelaxation + DiagonalMatrix: optimized path with + // lambdas + using PreconditionerType = DiagonalMatrix; + + using MyMatrixType = MySparseMatrix; + + MyMatrixType my_system_matrix(system_matrix); + + PreconditionRelaxation + preconditioner; + + PreconditionRelaxation::AdditionalData ad; + ad.relaxation = relaxation; + ad.n_iterations = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(my_system_matrix, ad); + + results.emplace_back(test(preconditioner, src, false)); + } + { // Test PreconditionRelaxation + DiagonalMatrix: optimized path using PreconditionerType = DiagonalMatrix;