From: Peter Munch Date: Wed, 14 Sep 2022 09:56:16 +0000 (+0200) Subject: PreconditionChebyshev: add specialization for A and P with ranges X-Git-Tag: v9.5.0-rc1~963^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F14271%2Fhead;p=dealii.git PreconditionChebyshev: add specialization for A and P with ranges --- diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 5998ecc035..d580e7d4b6 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -3046,14 +3046,21 @@ namespace internal // We need to have a separate declaration for static const members + // general case and the case that the preconditioner can work on + // ranges (covered by vector_updates()) template < typename MatrixType, typename VectorType, typename PreconditionerType, - std::enable_if_t, - int> * = nullptr> + std::enable_if_t< + !has_vmult_with_std_functions && + !(has_vmult_with_std_functions_for_precondition && + has_vmult_with_std_functions_for_precondition), + int> * = nullptr> inline void vmult_and_update(const MatrixType & matrix, const PreconditionerType &preconditioner, @@ -3079,6 +3086,126 @@ namespace internal solution); } + // case that both the operator and the preconditioner can work on + // subranges + template < + typename MatrixType, + typename VectorType, + typename PreconditionerType, + std::enable_if_t< + !has_vmult_with_std_functions && + (has_vmult_with_std_functions_for_precondition && + has_vmult_with_std_functions_for_precondition), + int> * = 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) + { + using Number = typename VectorType::value_type; + + const Number factor1 = factor1_; + const Number factor1_plus_1 = 1. + factor1_; + const Number factor2 = factor2_; + + if (iteration_index == 0) + { + preconditioner.vmult( + solution, + rhs, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'solution' before running the vmult operation + if (end_range > start_range) + std::memset(solution.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int start_range, const unsigned int end_range) { + const auto solution_ptr = solution.begin(); + + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + solution_ptr[i] *= factor2; + }); + } + else + { + temp_vector1.reinit(rhs, true); + temp_vector2.reinit(rhs, true); + + // 1) compute rediduum (including operator application) + 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) { + const auto rhs_ptr = rhs.begin(); + const auto tmp_ptr = temp_vector1.begin(); + + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + tmp_ptr[i] = factor2 * (rhs_ptr[i] - tmp_ptr[i]); + }); + + // 2) perform vector updates (including preconditioner application) + preconditioner.vmult( + temp_vector2, + temp_vector1, + [&](const unsigned int start_range, const unsigned int end_range) { + // zero 'temp_vector2' before running the vmult + // operation + if (end_range > start_range) + std::memset(temp_vector2.begin() + start_range, + 0, + sizeof(Number) * (end_range - start_range)); + }, + [&](const unsigned int start_range, const unsigned int end_range) { + const auto solution_ptr = solution.begin(); + const auto solution_old_ptr = solution_old.begin(); + const auto tmp_ptr = temp_vector2.begin(); + + if (iteration_index == 1) + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + tmp_ptr[i] = + factor1_plus_1 * solution_ptr[i] + factor2 * tmp_ptr[i]; + } + else + { + DEAL_II_OPENMP_SIMD_PRAGMA + for (std::size_t i = start_range; i < end_range; ++i) + tmp_ptr[i] = factor1_plus_1 * solution_ptr[i] - + factor1 * solution_old_ptr[i] + + factor2 * tmp_ptr[i]; + } + }); + + solution.swap(temp_vector2); + solution_old.swap(temp_vector2); + } + } + + // case that the operator can work on subranges and the preconditioner + // is a diagonal template +#include + +#include +#include + +#include +#include + +#include + +#include +#include +#include +#include + +#include + +#include "../tests.h" + +template +class MyDiagonalMatrix +{ +public: + void + vmult(VectorType &dst, const VectorType &src) const + { + diagonal_matrix.vmult(dst, src); + } + + void + Tvmult(VectorType &dst, const VectorType &src) const + { + diagonal_matrix.vmult(dst, src); + } + + VectorType & + get_vector() + { + return diagonal_matrix.get_vector(); + } + +private: + DiagonalMatrix diagonal_matrix; +}; + +template +class MyDiagonalMatrixWithPreAndPost +{ +public: + void + vmult(VectorType & dst, + const VectorType &src, + const std::function + &operation_before_matrix_vector_product = {}, + const std::function + &operation_after_matrix_vector_product = {}) const + { + if (operation_before_matrix_vector_product) + operation_before_matrix_vector_product(0, src.size()); + + diagonal_matrix.vmult(dst, src); + + if (operation_after_matrix_vector_product) + operation_after_matrix_vector_product(0, src.size()); + } + + VectorType & + get_vector() + { + return diagonal_matrix.get_vector(); + } + +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 + { + if (operation_before_matrix_vector_product) + operation_before_matrix_vector_product(0, src.size()); + + sparse_matrix.vmult(dst, src); + + if (operation_after_matrix_vector_product) + operation_after_matrix_vector_product(0, src.size()); + } + + types::global_dof_index + m() const + { + return sparse_matrix.m(); + } + + double + el(unsigned int i, unsigned int j) const + { + return sparse_matrix.el(i, j); + } + +private: + const SparseMatrixType &sparse_matrix; +}; + + +template +std::tuple +test(const PreconditionerType &preconditioner, const VectorType &src) +{ + VectorType dst; + dst.reinit(src); + + dst = 1.0; + preconditioner.vmult(dst, src); + const double norm_0 = dst.l2_norm(); + + dst = 1.0; + preconditioner.step(dst, src); + const double norm_1 = dst.l2_norm(); + + return std::tuple{norm_0, norm_1}; +} + + + +int +main() +{ + initlog(); + + using Number = double; + using VectorType = Vector; + using MatrixType = SparseMatrix; + + const unsigned int dim = 2; + const unsigned int degree = 1; + + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(3); + + FE_Q fe(degree); + QGauss quad(degree + 1); + MappingQ1 mapping; + + DoFHandler dof_handler(tria); + dof_handler.distribute_dofs(fe); + + MatrixType system_matrix; + + DynamicSparsityPattern dsp(dof_handler.n_dofs()); + DoFTools::make_sparsity_pattern(dof_handler, dsp); + + SparsityPattern sparsity_pattern; + sparsity_pattern.copy_from(dsp); + system_matrix.reinit(sparsity_pattern); + + MatrixCreator::create_laplace_matrix(mapping, + dof_handler, + quad, + system_matrix); + + VectorType diagonal(dof_handler.n_dofs()); + VectorType src(dof_handler.n_dofs()); + + for (const auto &entry : system_matrix) + if (entry.row() == entry.column()) + diagonal[entry.row()] = 1.0 / entry.value(); + + src = 1.0; + + std::vector n_iterationss{1, 2, 3}; + + for (const auto n_iterations : n_iterationss) + { + std::vector> results; + + { + // Test PreconditionChebyshev + DiagonalMatrix and matrix with + // pre/post + using PreconditionerType = DiagonalMatrix; + + using MyMatrixType = MySparseMatrix; + + MyMatrixType my_system_matrix(system_matrix); + + PreconditionChebyshev + preconditioner; + + PreconditionChebyshev:: + AdditionalData ad; + ad.degree = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(my_system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + /// Test PreconditionChebyshev + DiagonalMatrix and matrix without + // pre/post + using PreconditionerType = DiagonalMatrix; + + PreconditionChebyshev + preconditioner; + + PreconditionChebyshev:: + AdditionalData ad; + ad.degree = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + // Test PreconditionChebyshev + arbitrary preconditioner and matrix + // without pre/post + using PreconditionerType = MyDiagonalMatrix; + + PreconditionChebyshev + preconditioner; + + PreconditionChebyshev:: + AdditionalData ad; + ad.degree = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + // Test PreconditionChebyshev + preconditioner with pre/post and + // matrix without pre/post + using PreconditionerType = MyDiagonalMatrixWithPreAndPost; + + PreconditionChebyshev + preconditioner; + + PreconditionChebyshev:: + AdditionalData ad; + ad.degree = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + // Test PreconditionChebyshev + preconditioner with pre/post and + // matrix with pre/post + using PreconditionerType = MyDiagonalMatrixWithPreAndPost; + + using MyMatrixType = MySparseMatrix; + + MyMatrixType my_system_matrix(system_matrix); + + PreconditionChebyshev + preconditioner; + + PreconditionChebyshev:: + AdditionalData ad; + ad.degree = n_iterations; + ad.preconditioner = std::make_shared(); + ad.preconditioner->get_vector() = diagonal; + + preconditioner.initialize(my_system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + if (std::equal(results.begin(), + results.end(), + results.begin(), + [](const auto &a, const auto &b) { + if (std::abs(std::get<0>(a) - std::get<0>(b)) > 1e-6) + return false; + + if (std::abs(std::get<1>(a) - std::get<1>(b)) > 1e-6) + return false; + + return true; + })) + deallog << "OK!" << std::endl; + else + deallog << "ERROR!" << std::endl; + } +} diff --git a/tests/lac/precondition_chebyshev_07.output b/tests/lac/precondition_chebyshev_07.output new file mode 100644 index 0000000000..7ca52c43a9 --- /dev/null +++ b/tests/lac/precondition_chebyshev_07.output @@ -0,0 +1,4 @@ + +DEAL::OK! +DEAL::OK! +DEAL::OK!