#include <deal.II/base/utilities.h>
#include <deal.II/lac/affine_constraints.h>
+#include <deal.II/lac/block_vector_base.h>
#include <deal.II/lac/diagonal_matrix.h>
#include <deal.II/lac/identity_matrix.h>
#include <deal.II/lac/solver_cg.h>
/**
* 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.
*/
*/
double relaxation;
+ /**
+ * Number of smoothing steps to be performed.
+ */
+ unsigned int n_iterations;
+
/*
* Preconditioner.
*/
{
namespace PreconditionRelaxation
{
+ template <typename T, typename VectorType>
+ struct has_step
+ {
+ private:
+ static bool
+ detect(...);
+
+ template <typename U>
+ static decltype(
+ std::declval<U const>().step(std::declval<VectorType &>(),
+ std::declval<const VectorType &>()))
+ detect(const U &);
+
+ public:
+ static const bool value =
+ !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+ };
+
+ template <typename T, typename VectorType>
+ const bool has_step<T, VectorType>::value;
+
+ template <typename T, typename VectorType>
+ struct has_Tstep
+ {
+ private:
+ static bool
+ detect(...);
+
+ template <typename U>
+ static decltype(
+ std::declval<U const>().Tstep(std::declval<VectorType &>(),
+ std::declval<const VectorType &>()))
+ detect(const U &);
+
+ public:
+ static const bool value =
+ !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+ };
+
+ template <typename T, typename VectorType>
+ const bool has_Tstep<T, VectorType>::value;
+
+ template <typename T, typename VectorType>
+ struct has_jacobi_step
+ {
+ private:
+ static bool
+ detect(...);
+
+ template <typename U>
+ static decltype(
+ std::declval<U const>().Jacobi_step(std::declval<VectorType &>(),
+ std::declval<const VectorType &>(),
+ 1.0))
+ detect(const U &);
+
+ public:
+ static const bool value =
+ !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+ };
+
+ template <typename T, typename VectorType>
+ const bool has_jacobi_step<T, VectorType>::value;
+
+ template <typename T, typename VectorType>
+ struct has_SOR_step
+ {
+ private:
+ static bool
+ detect(...);
+
+ template <typename U>
+ static decltype(
+ std::declval<U const>().SOR_step(std::declval<VectorType &>(),
+ std::declval<const VectorType &>(),
+ 1.0))
+ detect(const U &);
+
+ public:
+ static const bool value =
+ !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+ };
+
+ template <typename T, typename VectorType>
+ const bool has_SOR_step<T, VectorType>::value;
+
+ template <typename T, typename VectorType>
+ struct has_SSOR_step
+ {
+ private:
+ static bool
+ detect(...);
+
+ template <typename U>
+ static decltype(
+ std::declval<U const>().SSOR_step(std::declval<VectorType &>(),
+ std::declval<const VectorType &>(),
+ 1.0))
+ detect(const U &);
+
+ public:
+ static const bool value =
+ !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
+ };
+
+ template <typename T, typename VectorType>
+ const bool has_SSOR_step<T, VectorType>::value;
+
template <typename MatrixType>
class PreconditionJacobiImpl
{
this->vmult(dst, src);
}
- template <typename VectorType>
+ template <
+ typename VectorType,
+ typename std::enable_if<has_jacobi_step<MatrixType, VectorType>::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<!has_jacobi_step<MatrixType, VectorType>::value,
+ MatrixType>::type * = nullptr>
+ void
+ step(VectorType &, const VectorType &) const
+ {
+ Assert(false,
+ ExcMessage(
+ "Matrix A does not provide a Jacobi_step() function!"));
+ }
+
template <typename VectorType>
void
Tstep(VectorType &dst, const VectorType &src) const
this->A->precondition_TSOR(dst, src, this->relaxation);
}
- template <typename VectorType>
+ template <
+ typename VectorType,
+ typename std::enable_if<has_SOR_step<MatrixType, VectorType>::value,
+ MatrixType>::type * = nullptr>
void
step(VectorType &dst, const VectorType &src) const
{
this->A->SOR_step(dst, src, this->relaxation);
}
- template <typename VectorType>
+ template <
+ typename VectorType,
+ typename std::enable_if<!has_SOR_step<MatrixType, VectorType>::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<has_SOR_step<MatrixType, VectorType>::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<!has_SOR_step<MatrixType, VectorType>::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<const MatrixType> A;
const double relaxation;
pos_right_of_diagonal);
}
- template <typename VectorType>
+ template <
+ typename VectorType,
+ typename std::enable_if<has_SSOR_step<MatrixType, VectorType>::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<!has_SSOR_step<MatrixType, VectorType>::value,
+ MatrixType>::type * = nullptr>
+ void
+ step(VectorType &, const VectorType &) const
+ {
+ Assert(false,
+ ExcMessage("Matrix A does not provide a SSOR_step() function!"));
+ }
+
template <typename VectorType>
void
Tstep(VectorType &dst, const VectorType &src) const
const std::vector<size_type> &inverse_permutation;
};
- template <typename T, typename VectorType>
- struct has_step
- {
- private:
- static bool
- detect(...);
-
- template <typename U>
- static decltype(
- std::declval<U const>().step(std::declval<VectorType &>(),
- std::declval<const VectorType &>()))
- detect(const U &);
-
- public:
- static const bool value =
- !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
- };
-
- template <typename T, typename VectorType>
- const bool has_step<T, VectorType>::value;
-
- template <typename T, typename VectorType>
- struct has_Tstep
- {
- private:
- static bool
- detect(...);
-
- template <typename U>
- static decltype(
- std::declval<U const>().Tstep(std::declval<VectorType &>(),
- std::declval<const VectorType &>()))
- detect(const U &);
-
- public:
- static const bool value =
- !std::is_same<bool, decltype(detect(std::declval<T>()))>::value;
- };
-
- template <typename T, typename VectorType>
- const bool has_Tstep<T, VectorType>::value;
-
template <
+ typename MatrixType,
typename PreconditionerType,
typename VectorType,
typename std::enable_if<has_step<PreconditionerType, VectorType>::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<!has_step<PreconditionerType, VectorType>::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<has_Tstep<PreconditionerType, VectorType>::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<!has_Tstep<PreconditionerType, VectorType>::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 <typename MatrixType,
+ typename PreconditionerType,
+ typename VectorType>
+ 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 <typename MatrixType,
+ typename VectorType,
+ typename std::enable_if<!IsBlockVector<VectorType>::value,
+ VectorType>::type * = nullptr>
+ void
+ step_operations(const MatrixType & A,
+ const DiagonalMatrix<VectorType> &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
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;
}
inline void
PreconditionRelaxation<MatrixType, PreconditionerType>::clear()
{
- A = nullptr;
+ A = nullptr;
+ preconditioner = nullptr;
}
template <typename MatrixType, typename PreconditionerType>
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 <typename MatrixType, typename PreconditionerType>
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 <typename MatrixType, typename PreconditionerType>
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 <typename MatrixType, typename PreconditionerType>
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);
}
//---------------------------------------------------------------------------
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<PreconditionerType>(A, parameters_in.relaxation);
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<PreconditionerType>(A, parameters_in.relaxation);
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<PreconditionerType>(A, parameters_in.relaxation);
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<PreconditionerType>(A, parameters_in.relaxation, p, ip);
template <typename MatrixType, typename PreconditionerType>
inline PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData::
- AdditionalData(const double relaxation)
+ AdditionalData(const double relaxation, const unsigned int n_iterations)
: relaxation(relaxation)
+ , n_iterations(n_iterations)
{}
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <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 PreconditionerType, typename VectorType>
+std::tuple<double, double, double, double>
+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<double, double, double, double>{norm_0,
+ norm_1,
+ norm_2,
+ norm_3};
+}
+
+
+
+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<double> relaxations{1.0, 0.9, 1.1};
+ std::vector<unsigned int> n_iterationss{1, 2, 3};
+
+ for (const auto relaxation : relaxations)
+ for (const auto n_iterations : n_iterationss)
+ {
+ std::vector<std::tuple<double, double, double, double>> results;
+
+ {
+ // Test PreconditionJacobi
+ PreconditionJacobi<MatrixType> preconditioner;
+
+ PreconditionJacobi<MatrixType>::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<VectorType>;
+
+ PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
+
+ PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData
+ ad;
+ ad.relaxation = relaxation;
+ ad.n_iterations = 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 PreconditionRelaxation + wrapper around DiagonalMatrix:
+ // optimized path cannot be used
+ using PreconditionerType = MyDiagonalMatrix<VectorType>;
+
+ PreconditionRelaxation<MatrixType, PreconditionerType> preconditioner;
+
+ PreconditionRelaxation<MatrixType, PreconditionerType>::AdditionalData
+ ad;
+ ad.relaxation = relaxation;
+ ad.n_iterations = n_iterations;
+ ad.preconditioner = std::make_shared<PreconditionerType>();
+ 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;
+ }
+}
projects / dealii.git / commitdiff