From 8b62d1ff5f5015c24d8302ad4dbbd68162fd902e Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Fri, 31 Dec 2021 16:30:21 +0100 Subject: [PATCH] Extand PreconditionRelaxation --- doc/news/changes/minor/20220101Munch | 8 + include/deal.II/lac/precondition.h | 439 ++++++++++++++++---- include/deal.II/matrix_free/operators.h | 1 - tests/lac/precondition_relaxation_01.cc | 221 ++++++++++ tests/lac/precondition_relaxation_01.output | 10 + 5 files changed, 599 insertions(+), 80 deletions(-) create mode 100644 doc/news/changes/minor/20220101Munch create mode 100644 tests/lac/precondition_relaxation_01.cc create mode 100644 tests/lac/precondition_relaxation_01.output diff --git a/doc/news/changes/minor/20220101Munch b/doc/news/changes/minor/20220101Munch new file mode 100644 index 0000000000..a10abb7e98 --- /dev/null +++ b/doc/news/changes/minor/20220101Munch @@ -0,0 +1,8 @@ +Improved: The class PreconditionRelaxation has been refactored. +Similarly to PreconditionChebyshev, it can be now used +stand-alone and users can attach their own preconditioners. +Furthermore, users can specify multiple iteration steps, +which is particularly useful if PreconditionRelaxation is +used as smoother in a multigrid algorithm. +
+(Peter Munch, 2022/01/01) diff --git a/include/deal.II/lac/precondition.h b/include/deal.II/lac/precondition.h index 4b0377e8a3..cc2a1ac3fb 100644 --- a/include/deal.II/lac/precondition.h +++ b/include/deal.II/lac/precondition.h @@ -29,6 +29,7 @@ #include #include +#include #include #include #include @@ -418,13 +419,20 @@ public: /** * Constructor. */ - AdditionalData(const double relaxation = 1.); + AdditionalData(const double relaxation = 1., + const unsigned int n_iterations = 1); /** * Relaxation parameter. */ double relaxation; + /** + * Number of smoothing steps to be performed. + */ + unsigned int n_iterations; + + /* * Preconditioner. */ @@ -500,6 +508,11 @@ protected: */ double relaxation; + /** + * Number of smoothing steps to be performed. + */ + unsigned int n_iterations; + /* * Preconditioner. */ @@ -514,6 +527,114 @@ namespace internal { namespace PreconditionRelaxation { + template + struct has_step + { + private: + static bool + detect(...); + + template + static decltype( + std::declval().step(std::declval(), + std::declval())) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + template + const bool has_step::value; + + template + struct has_Tstep + { + private: + static bool + detect(...); + + template + static decltype( + std::declval().Tstep(std::declval(), + std::declval())) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + template + const bool has_Tstep::value; + + template + struct has_jacobi_step + { + private: + static bool + detect(...); + + template + static decltype( + std::declval().Jacobi_step(std::declval(), + std::declval(), + 1.0)) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + template + const bool has_jacobi_step::value; + + template + struct has_SOR_step + { + private: + static bool + detect(...); + + template + static decltype( + std::declval().SOR_step(std::declval(), + std::declval(), + 1.0)) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + template + const bool has_SOR_step::value; + + template + struct has_SSOR_step + { + private: + static bool + detect(...); + + template + static decltype( + std::declval().SSOR_step(std::declval(), + std::declval(), + 1.0)) + detect(const U &); + + public: + static const bool value = + !std::is_same()))>::value; + }; + + template + const bool has_SSOR_step::value; + template class PreconditionJacobiImpl { @@ -538,13 +659,28 @@ namespace internal this->vmult(dst, src); } - template + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> void step(VectorType &dst, const VectorType &src) const { this->A->Jacobi_step(dst, src, this->relaxation); } + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> + void + step(VectorType &, const VectorType &) const + { + Assert(false, + ExcMessage( + "Matrix A does not provide a Jacobi_step() function!")); + } + template void Tstep(VectorType &dst, const VectorType &src) const @@ -581,20 +717,48 @@ namespace internal this->A->precondition_TSOR(dst, src, this->relaxation); } - template + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> void step(VectorType &dst, const VectorType &src) const { this->A->SOR_step(dst, src, this->relaxation); } - template + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> + void + step(VectorType &, const VectorType &) const + { + Assert(false, + ExcMessage("Matrix A does not provide a SOR_step() function!")); + } + + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> void Tstep(VectorType &dst, const VectorType &src) const { this->A->TSOR_step(dst, src, this->relaxation); } + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> + void + Tstep(VectorType &, const VectorType &) const + { + Assert(false, + ExcMessage("Matrix A does not provide a TSOR_step() function!")); + } + private: const SmartPointer A; const double relaxation; @@ -658,13 +822,27 @@ namespace internal pos_right_of_diagonal); } - template + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> void step(VectorType &dst, const VectorType &src) const { this->A->SSOR_step(dst, src, this->relaxation); } + template < + typename VectorType, + typename std::enable_if::value, + MatrixType>::type * = nullptr> + void + step(VectorType &, const VectorType &) const + { + Assert(false, + ExcMessage("Matrix A does not provide a SSOR_step() function!")); + } + template void Tstep(VectorType &dst, const VectorType &src) const @@ -724,104 +902,175 @@ namespace internal const std::vector &inverse_permutation; }; - template - struct has_step - { - private: - static bool - detect(...); - - template - static decltype( - std::declval().step(std::declval(), - std::declval())) - detect(const U &); - - public: - static const bool value = - !std::is_same()))>::value; - }; - - template - const bool has_step::value; - - template - struct has_Tstep - { - private: - static bool - detect(...); - - template - static decltype( - std::declval().Tstep(std::declval(), - std::declval())) - detect(const U &); - - public: - static const bool value = - !std::is_same()))>::value; - }; - - template - const bool has_Tstep::value; - template < + typename MatrixType, typename PreconditionerType, typename VectorType, typename std::enable_if::value, PreconditionerType>::type * = nullptr> void - step(const PreconditionerType &preconditioner, + step(const MatrixType &, + const PreconditionerType &preconditioner, VectorType & dst, - const VectorType & src) + const VectorType & src, + const double relaxation, + VectorType &, + VectorType &) { + Assert(relaxation == 1.0, ExcInternalError()); + + (void)relaxation; + preconditioner.step(dst, src); } template < + typename MatrixType, typename PreconditionerType, typename VectorType, typename std::enable_if::value, PreconditionerType>::type * = nullptr> void - step(const PreconditionerType &preconditioner, + step(const MatrixType & A, + const PreconditionerType &preconditioner, VectorType & dst, - const VectorType & src) + const VectorType & src, + const double relaxation, + VectorType & residual, + VectorType & tmp) { - Assert(false, ExcNotImplemented()); - (void)preconditioner; - (void)dst; - (void)src; + residual.reinit(dst, true); + tmp.reinit(dst, true); + + A.vmult(residual, dst); + residual.sadd(-1.0, 1.0, src); + + preconditioner.vmult(tmp, residual); + dst.add(relaxation, tmp); } template < + typename MatrixType, typename PreconditionerType, typename VectorType, typename std::enable_if::value, PreconditionerType>::type * = nullptr> void - Tstep(const PreconditionerType &preconditioner, + Tstep(const MatrixType &, + const PreconditionerType &preconditioner, VectorType & dst, - const VectorType & src) + const VectorType & src, + const double relaxation, + VectorType &, + VectorType &) { + Assert(relaxation == 1.0, ExcInternalError()); + + (void)relaxation; + preconditioner.Tstep(dst, src); } template < + typename MatrixType, typename PreconditionerType, typename VectorType, typename std::enable_if::value, PreconditionerType>::type * = nullptr> void - Tstep(const PreconditionerType &preconditioner, + Tstep(const MatrixType & A, + const PreconditionerType &preconditioner, VectorType & dst, - const VectorType & src) + const VectorType & src, + const double relaxation, + VectorType & residual, + VectorType & tmp) { - Assert(false, ExcNotImplemented()); - (void)preconditioner; - (void)dst; - (void)src; + residual.reinit(dst, true); + tmp.reinit(dst, true); + + A.Tvmult(residual, dst); + residual.sadd(-1.0, 1.0, src); + + preconditioner.Tvmult(tmp, residual); + dst.add(relaxation, tmp); + } + + template + void + step_operations(const MatrixType & A, + const PreconditionerType &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & tmp1, + VectorType & tmp2, + const unsigned int i, + const bool transposed) + { + if (i == 0) + { + if (transposed) + preconditioner.Tvmult(dst, src); + else + preconditioner.vmult(dst, src); + + if (relaxation != 1.0) + dst *= relaxation; + } + else + { + if (transposed) + Tstep(A, preconditioner, dst, src, relaxation, tmp1, tmp2); + else + step(A, preconditioner, dst, src, relaxation, tmp1, tmp2); + } + } + + template ::value, + VectorType>::type * = nullptr> + void + step_operations(const MatrixType & A, + const DiagonalMatrix &preconditioner, + VectorType & dst, + const VectorType & src, + const double relaxation, + VectorType & residual, + VectorType &, + const unsigned int i, + const bool transposed) + { + if (i == 0) + { + const auto dst_ptr = dst.begin(); + const auto src_ptr = src.begin(); + const auto 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]; + } + else + { + residual.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(); + + if (transposed) + A.Tvmult(residual, dst); + else + A.vmult(residual, 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]; + } } } // namespace PreconditionRelaxation @@ -1752,9 +2001,13 @@ PreconditionRelaxation::initialize( const MatrixType & rA, const AdditionalData ¶meters) { - A = &rA; - relaxation = parameters.relaxation; + A = &rA; + relaxation = parameters.relaxation; + + Assert(parameters.preconditioner, ExcNotInitialized()); + preconditioner = parameters.preconditioner; + n_iterations = parameters.n_iterations; } @@ -1762,7 +2015,8 @@ template inline void PreconditionRelaxation::clear() { - A = nullptr; + A = nullptr; + preconditioner = nullptr; } template @@ -1790,10 +2044,14 @@ PreconditionRelaxation::vmult( VectorType & dst, const VectorType &src) const { - preconditioner->vmult(dst, src); + Assert(this->A != nullptr, ExcNotInitialized()); + Assert(this->preconditioner != nullptr, ExcNotInitialized()); - if (this->relaxation != 1.0) - dst *= this->relaxation; + VectorType tmp1, tmp2; + + for (unsigned int i = 0; i < n_iterations; ++i) + internal::PreconditionRelaxation::step_operations( + *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false); } template @@ -1803,10 +2061,14 @@ PreconditionRelaxation::Tvmult( VectorType & dst, const VectorType &src) const { - preconditioner->Tvmult(dst, src); + Assert(this->A != nullptr, ExcNotInitialized()); + Assert(this->preconditioner != nullptr, ExcNotInitialized()); + + VectorType tmp1, tmp2; - if (this->relaxation != 1.0) - dst *= this->relaxation; + for (unsigned int i = 0; i < n_iterations; ++i) + internal::PreconditionRelaxation::step_operations( + *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true); } template @@ -1816,7 +2078,14 @@ PreconditionRelaxation::step( VectorType & dst, const VectorType &src) const { - internal::PreconditionRelaxation::step(*preconditioner, dst, src); + Assert(this->A != nullptr, ExcNotInitialized()); + Assert(this->preconditioner != nullptr, ExcNotInitialized()); + + VectorType tmp1, tmp2; + + for (unsigned int i = 1; i <= n_iterations; ++i) + internal::PreconditionRelaxation::step_operations( + *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, false); } template @@ -1826,7 +2095,14 @@ PreconditionRelaxation::Tstep( VectorType & dst, const VectorType &src) const { - internal::PreconditionRelaxation::Tstep(*preconditioner, dst, src); + Assert(this->A != nullptr, ExcNotInitialized()); + Assert(this->preconditioner != nullptr, ExcNotInitialized()); + + VectorType tmp1, tmp2; + + for (unsigned int i = 1; i <= n_iterations; ++i) + internal::PreconditionRelaxation::step_operations( + *A, *preconditioner, dst, src, relaxation, tmp1, tmp2, i, true); } //--------------------------------------------------------------------------- @@ -1839,7 +2115,8 @@ PreconditionJacobi::initialize(const MatrixType & A, Assert(parameters_in.preconditioner == nullptr, ExcInternalError()); AdditionalData parameters; - parameters.relaxation = 1.0; + parameters.relaxation = 1.0; + parameters.n_iterations = parameters_in.n_iterations; parameters.preconditioner = std::make_shared(A, parameters_in.relaxation); @@ -1856,7 +2133,8 @@ PreconditionSOR::initialize(const MatrixType & A, Assert(parameters_in.preconditioner == nullptr, ExcInternalError()); AdditionalData parameters; - parameters.relaxation = 1.0; + parameters.relaxation = 1.0; + parameters.n_iterations = parameters_in.n_iterations; parameters.preconditioner = std::make_shared(A, parameters_in.relaxation); @@ -1873,7 +2151,8 @@ PreconditionSSOR::initialize(const MatrixType & A, Assert(parameters_in.preconditioner == nullptr, ExcInternalError()); AdditionalData parameters; - parameters.relaxation = 1.0; + parameters.relaxation = 1.0; + parameters.n_iterations = parameters_in.n_iterations; parameters.preconditioner = std::make_shared(A, parameters_in.relaxation); @@ -1895,7 +2174,8 @@ PreconditionPSOR::initialize( Assert(parameters_in.preconditioner == nullptr, ExcInternalError()); typename BaseClass::AdditionalData parameters; - parameters.relaxation = 1.0; + parameters.relaxation = 1.0; + parameters.n_iterations = parameters_in.n_iterations; parameters.preconditioner = std::make_shared(A, parameters_in.relaxation, p, ip); @@ -1952,8 +2232,9 @@ PreconditionUseMatrix::vmult( template inline PreconditionRelaxation::AdditionalData:: - AdditionalData(const double relaxation) + AdditionalData(const double relaxation, const unsigned int n_iterations) : relaxation(relaxation) + , n_iterations(n_iterations) {} diff --git a/include/deal.II/matrix_free/operators.h b/include/deal.II/matrix_free/operators.h index 64f18fae96..dd618bcd9a 100644 --- a/include/deal.II/matrix_free/operators.h +++ b/include/deal.II/matrix_free/operators.h @@ -385,7 +385,6 @@ namespace MatrixFreeOperators const std::shared_ptr> & get_matrix_diagonal() const; - /** * Apply the Jacobi preconditioner, which multiplies every element of the * src vector by the inverse of the respective diagonal element and diff --git a/tests/lac/precondition_relaxation_01.cc b/tests/lac/precondition_relaxation_01.cc new file mode 100644 index 0000000000..cc0a89061c --- /dev/null +++ b/tests/lac/precondition_relaxation_01.cc @@ -0,0 +1,221 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2021 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test PreconditionRelaxation for different Jacobi preconditioners. + + +#include +#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 +std::tuple +test(const PreconditionerType &preconditioner, const VectorType &src) +{ + VectorType dst; + dst.reinit(src); + + preconditioner.vmult(dst, src); + const double norm_0 = dst.l2_norm(); + + 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); + const double norm_3 = dst.l2_norm(); + + return std::tuple{norm_0, + norm_1, + norm_2, + norm_3}; +} + + + +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 relaxations{1.0, 0.9, 1.1}; + std::vector n_iterationss{1, 2, 3}; + + for (const auto relaxation : relaxations) + for (const auto n_iterations : n_iterationss) + { + std::vector> results; + + { + // Test PreconditionJacobi + PreconditionJacobi preconditioner; + + PreconditionJacobi::AdditionalData ad; + ad.relaxation = relaxation; + ad.n_iterations = n_iterations; + + preconditioner.initialize(system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + // Test PreconditionRelaxation + DiagonalMatrix: optimized path + using PreconditionerType = DiagonalMatrix; + + 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(system_matrix, ad); + + results.emplace_back(test(preconditioner, src)); + } + + { + // Test PreconditionRelaxation + wrapper around DiagonalMatrix: + // optimized path cannot be used + using PreconditionerType = MyDiagonalMatrix; + + 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(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; + + if (std::abs(std::get<2>(a) - std::get<2>(b)) > 1e-6) + return false; + + if (std::abs(std::get<3>(a) - std::get<3>(b)) > 1e-6) + return false; + + return true; + })) + deallog << "OK!" << std::endl; + else + deallog << "ERROR!" << std::endl; + } +} diff --git a/tests/lac/precondition_relaxation_01.output b/tests/lac/precondition_relaxation_01.output new file mode 100644 index 0000000000..f738e07d35 --- /dev/null +++ b/tests/lac/precondition_relaxation_01.output @@ -0,0 +1,10 @@ + +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! +DEAL::OK! -- 2.39.5