// 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<!has_vmult_with_std_functions<MatrixType,
- VectorType,
- PreconditionerType>,
- int> * = nullptr>
+ std::enable_if_t<
+ !has_vmult_with_std_functions<MatrixType,
+ VectorType,
+ PreconditionerType> &&
+ !(has_vmult_with_std_functions_for_precondition<PreconditionerType,
+ VectorType> &&
+ has_vmult_with_std_functions_for_precondition<MatrixType,
+ VectorType>),
+ int> * = nullptr>
inline void
vmult_and_update(const MatrixType & matrix,
const PreconditionerType &preconditioner,
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<MatrixType,
+ VectorType,
+ PreconditionerType> &&
+ (has_vmult_with_std_functions_for_precondition<PreconditionerType,
+ VectorType> &&
+ has_vmult_with_std_functions_for_precondition<MatrixType,
+ VectorType>),
+ 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 <typename MatrixType,
typename VectorType,
typename PreconditionerType,
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2021 - 2022 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 PreconditionChebyshev for different optimization levels. The
+// test is similar to precondition_relaxation_01.
+
+
+#include <deal.II/base/mpi.h>
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q1.h>
+
+#include <deal.II/grid/grid_generator.h>
+
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/lac/vector.h>
+
+#include <deal.II/numerics/matrix_tools.h>
+
+#include "../tests.h"
+
+template <typename VectorType>
+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<VectorType> diagonal_matrix;
+};
+
+template <typename VectorType>
+class MyDiagonalMatrixWithPreAndPost
+{
+public:
+ void
+ vmult(VectorType & dst,
+ const VectorType &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_matrix_vector_product = {},
+ const std::function<void(const unsigned int, const unsigned int)>
+ &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<VectorType> diagonal_matrix;
+};
+
+template <typename SparseMatrixType>
+class MySparseMatrix : public Subscriptor
+{
+public:
+ MySparseMatrix(const SparseMatrixType &sparse_matrix)
+ : sparse_matrix(sparse_matrix)
+ {}
+
+ template <typename VectorType>
+ void
+ vmult(VectorType & dst,
+ const VectorType &src,
+ const std::function<void(const unsigned int, const unsigned int)>
+ &operation_before_matrix_vector_product = {},
+ const std::function<void(const unsigned int, const unsigned int)>
+ &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 <typename PreconditionerType, typename VectorType>
+std::tuple<double, double>
+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<double, double>{norm_0, norm_1};
+}
+
+
+
+int
+main()
+{
+ initlog();
+
+ using Number = double;
+ using VectorType = Vector<Number>;
+ using MatrixType = SparseMatrix<Number>;
+
+ const unsigned int dim = 2;
+ const unsigned int degree = 1;
+
+ Triangulation<dim> tria;
+ GridGenerator::hyper_cube(tria);
+ tria.refine_global(3);
+
+ FE_Q<dim> fe(degree);
+ QGauss<dim> quad(degree + 1);
+ MappingQ1<dim> mapping;
+
+ DoFHandler<dim> 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<unsigned int> n_iterationss{1, 2, 3};
+
+ for (const auto n_iterations : n_iterationss)
+ {
+ std::vector<std::tuple<double, double>> results;
+
+ {
+ // Test PreconditionChebyshev + DiagonalMatrix and matrix with
+ // pre/post
+ using PreconditionerType = DiagonalMatrix<VectorType>;
+
+ using MyMatrixType = MySparseMatrix<MatrixType>;
+
+ MyMatrixType my_system_matrix(system_matrix);
+
+ PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>
+ preconditioner;
+
+ PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>::
+ AdditionalData ad;
+ ad.degree = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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<VectorType>;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+ preconditioner;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+ AdditionalData ad;
+ ad.degree = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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<VectorType>;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+ preconditioner;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+ AdditionalData ad;
+ ad.degree = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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<VectorType>;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>
+ preconditioner;
+
+ PreconditionChebyshev<MatrixType, VectorType, PreconditionerType>::
+ AdditionalData ad;
+ ad.degree = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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<VectorType>;
+
+ using MyMatrixType = MySparseMatrix<MatrixType>;
+
+ MyMatrixType my_system_matrix(system_matrix);
+
+ PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>
+ preconditioner;
+
+ PreconditionChebyshev<MyMatrixType, VectorType, PreconditionerType>::
+ AdditionalData ad;
+ ad.degree = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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;
+ }
+}
projects / dealii.git / commitdiff