From: Matthias Maier Date: Thu, 9 Apr 2015 00:40:54 +0000 (+0200) Subject: Generalize the "linop" constructor X-Git-Tag: v8.3.0-rc1~242^2~13 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=8febacf17dfec6a487d2456e245f937a959e710b;p=dealii.git Generalize the "linop" constructor --- diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index d95e21dcf2..e437d30402 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -33,10 +33,17 @@ template class Vector; template class BlockVector; template class LinearOperator; + +template , + typename Domain = Range, + typename Exemplar, + typename Matrix> +LinearOperator linop (const Exemplar &, const Matrix &); + template , typename Domain = Range, typename Matrix> -LinearOperator linop(const Matrix &matrix); +LinearOperator linop (const Matrix &); /** * A class to store the abstract concept of a linear operator. @@ -744,6 +751,120 @@ public: }; +namespace +{ + // A trait class that determines whether type T provides public + // (templated or non-templated) vmult_add and Tvmult_add member functions + + template + class has_vmult_add + { + template + static std::false_type test(...); + + template + static std::true_type test(decltype(&C::vmult_add), + decltype(&C::Tvmult_add)); + + template + static std::true_type test(decltype(&C::template vmult_add), + decltype(&C::template Tvmult_add)); + + template + static std::true_type test(decltype(&C::template vmult_add), + decltype(&C::template Tvmult_add)); + + public: + // type is std::true_type if Matrix provides vmult_add and Tvmult_add, + // otherwise it is std::false_type + + typedef decltype(test(0, 0)) type; + }; + + + // A helper class to add the full matrix interface to a LinearOperator + + template + class MatrixInterfaceWithVmultAdd + { + public: + + template + void operator()(LinearOperator &op, const Matrix &matrix) + { + op.vmult = [&matrix](Range &v, const Domain &u) + { + matrix.vmult(v,u); + }; + + op.vmult_add = [&matrix](Range &v, const Domain &u) + { + matrix.vmult_add(v,u); + }; + + op.Tvmult = [&matrix](Domain &v, const Range &u) + { + matrix.Tvmult(v,u); + }; + + op.Tvmult_add = [&matrix](Domain &v, const Range &u) + { + matrix.Tvmult_add(v,u); + }; + } + }; + + + // A helper class to add a reduced matrix interface to a LinearOperator + // (typically provided by Preconditioner classes) + + template + class MatrixInterfaceWithoutVmultAdd + { + public: + + template + void operator()(LinearOperator &op, const Matrix &matrix) + { + op.vmult = [&matrix](Range &v, const Domain &u) + { + matrix.vmult(v,u); + }; + + op.vmult_add = [op](Range &v, const Domain &u) + { + // Capture op by value to have a valid GrowingVectorMemory object + // that is bound to the lifetime of this function object + + Range *i = op.range_vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/true); + op.vmult(*i, u); + v += *i; + op.range_vector_memory.free(i); + }; + + op.Tvmult = [&matrix](Domain &v, const Range &u) + { + matrix.Tvmult(v,u); + }; + + op.Tvmult_add = [op](Domain &v, const Range &u) + { + // Capture op by value to have a valid GrowingVectorMemory object + // that is bound to the lifetime of this function object + + Domain *i = op.domain_vector_memory.alloc(); + op.reinit_domain_vector(*i, /*bool fast =*/true); + op.Tvmult(*i, u); + v += *i; + op.domain_vector_memory.free(i); + }; + } + }; + +} /* namespace */ + + /** * A function that encapsulates generic @p matrix objects that act on a * compatible Vector type into a LinearOperator. The LinearOperator object @@ -769,60 +890,82 @@ public: * // Application of matrix to vector src, writes the result into dst. * vmult(VECTOR &dst, const VECTOR &src); * - * // Application of matrix to vector src, adds the result to dst. - * vmult_add(VECTOR &dst, const VECTOR &src); - * * // Application of the transpose of matrix to vector src, writes the * // result into dst. (Depending on the usage of the linear operator * // class this can be a dummy implementation throwing an error.) * Tvmult(VECTOR &dst, const VECTOR &src); + * }; + * @endcode + * + * The following (optional) interface is used if available: + * + * @code + * class Matrix + * { + * public: + * // Application of matrix to vector src, adds the result to dst. + * vmult_add(VECTOR &dst, const VECTOR &src); * * // Application of the transpose of matrix to vector src, adds the - * // result to dst. (Depending on the usage of the linear operator class - * // this can be a dummy implementation throwing an error.) + * // result to dst. * Tvmult_add(VECTOR &dst, const VECTOR &src); * }; * @endcode * + * If the Matrix does not provide vmult_add and + * Tvmult_add, they are implemented in terms of + * vmult and Tvmult (requiring intermediate + * storage). + * * @author: Matthias Maier, 2015 */ template , typename Domain = Range, typename Matrix> LinearOperator linop(const Matrix &matrix) +{ + // implement with the more generic variant below... + return linop(matrix, matrix); +} + + +/** + * Variant of above function that takes an Exemplar object @p exemplar as + * an additional reference. This object is used to populate the + * reinit_domain_vector and reinit_range_vector function objects. The + * reference @p matrix is used to construct vmult, Tvmult, etc. + * + * This variant can, for example, be used to encapsulate preconditioners + * (that typically do not expose any information about the underlying + * matrix). + * + * @author: Matthias Maier, 2015 + */ +template , + typename Domain = Range, + typename Exemplar, + typename Matrix> +LinearOperator linop(const Exemplar &exemplar, + const Matrix &matrix) { LinearOperator return_op; - // Always store a reference to matrix in the lambda functions. - // This ensures that a modification of the matrix after the creation of a - // LinearOperator wrapper is respected - further a matrix cannot usually - // be copied... + // Always store a reference to matrix and exemplar in the lambda + // functions. This ensures that a modification of the matrix after the + // creation of a LinearOperator wrapper is respected - further a matrix + // or an exemplar cannot usually be copied... return_op.reinit_range_vector = - ReinitRangeFactory().operator()(matrix); + ReinitRangeFactory().operator()(exemplar); return_op.reinit_domain_vector = - ReinitDomainFactory().operator()(matrix); - - return_op.vmult = [&matrix](Range &v, const Domain &u) - { - matrix.vmult(v,u); - }; - - return_op.vmult_add = [&matrix](Range &v, const Domain &u) - { - matrix.vmult_add(v,u); - }; + ReinitDomainFactory().operator()(exemplar); - return_op.Tvmult = [&matrix](Domain &v, const Range &u) - { - matrix.Tvmult(v,u); - }; - - return_op.Tvmult_add = [&matrix](Domain &v, const Range &u) - { - matrix.Tvmult_add(v,u); - }; + typename std::conditional< + has_vmult_add::type::value, + MatrixInterfaceWithVmultAdd, + MatrixInterfaceWithoutVmultAdd>::type(). + operator()(return_op, matrix); return return_op; } diff --git a/tests/lac/linear_operator_06.cc b/tests/lac/linear_operator_06.cc new file mode 100644 index 0000000000..393b548c61 --- /dev/null +++ b/tests/lac/linear_operator_06.cc @@ -0,0 +1,239 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2015 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + +// Test linop constructor for full and reduced interface: + +#include "../tests.h" + +#include +#include +#include + +class MyMatrix1 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + void vmult(Vector &, const Vector &) const + { + deallog << "MyMatrix1::vmult" << std::endl; + } + + void Tvmult(Vector &, const Vector &) const + { + deallog << "MyMatrix1::Tvmult" << std::endl; + } + +}; + +class MyMatrix2 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + void vmult(Vector &, const Vector &) const + { + deallog << "MyMatrix2::vmult" << std::endl; + } + + void Tvmult(Vector &, const Vector &) const + { + deallog << "MyMatrix2::Tvmult" << std::endl; + } + + void vmult_add(Vector &, const Vector &) const + { + deallog << "MyMatrix2::vmult_add" << std::endl; + } + + void Tvmult_add(Vector &, const Vector &) const + { + deallog << "MyMatrix2::Tvmult_add" << std::endl; + } +}; + +class MyMatrix3 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + template + void vmult(OutVector &, const InVector &) const + { + deallog << "MyMatrix3::vmult" << std::endl; + } + + template + void Tvmult(OutVector &, const InVector &) const + { + deallog << "MyMatrix3::Tvmult" << std::endl; + } + + template + void vmult_add(OutVector &, const InVector &) const + { + deallog << "MyMatrix3::vmult_add" << std::endl; + } + + template + void Tvmult_add(OutVector &, const InVector &) const + { + deallog << "MyMatrix3::Tvmult_add" << std::endl; + } +}; + +class MyMatrix4 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + template + void vmult(OutVector &, const InVector &, bool = true) const + { + deallog << "MyMatrix4::vmult" << std::endl; + } + + template + void Tvmult(OutVector &, const InVector &, bool = true) const + { + deallog << "MyMatrix4::Tvmult" << std::endl; + } + + template + void vmult_add(OutVector &, const InVector &, bool = true) const + { + deallog << "MyMatrix4::vmult_add" << std::endl; + } + + template + void Tvmult_add(OutVector &, const InVector &, bool = true) const + { + deallog << "MyMatrix4::Tvmult_add" << std::endl; + } +}; + +class MyMatrix5 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + template + void vmult(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix5::vmult" << std::endl; + } + + template + void Tvmult(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix5::Tvmult" << std::endl; + } + + template + void vmult_add(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix5::vmult_add" << std::endl; + } + + template + void Tvmult_add(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix5::Tvmult_add" << std::endl; + } +}; + +class MyMatrix6 +{ +public: + size_t m() const { return 1; }; + size_t n() const { return 1; }; + + template + void vmult(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix6::vmult" << std::endl; + } + + template + void Tvmult(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix6::Tvmult" << std::endl; + } + + template + void vmult_add(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix6::vmult_add" << std::endl; + } + + template + void Tvmult_add(Vector &, const Vector &, bool = true) const + { + deallog << "MyMatrix6::Tvmult_add" << std::endl; + } +}; + +using namespace dealii; + +int main() +{ + initlog(); + deallog << std::setprecision(10); + + Vector u (1); + Vector v (1); + + MyMatrix1 m1; + MyMatrix2 m2; + MyMatrix3 m3; + MyMatrix4 m4; + MyMatrix5 m5; + MyMatrix6 m6; + + linop(m1).vmult(v, u); + linop(m1).Tvmult(v, u); + linop(m1).vmult_add(v, u); + linop(m1).Tvmult_add(v, u); + + linop(m2).vmult(v, u); + linop(m2).Tvmult(v, u); + linop(m2).vmult_add(v, u); + linop(m2).Tvmult_add(v, u); + + linop(m3).vmult(v, u); + linop(m3).Tvmult(v, u); + linop(m3).vmult_add(v, u); + linop(m3).Tvmult_add(v, u); + + linop(m4).vmult(v, u); + linop(m4).Tvmult(v, u); + linop(m4).vmult_add(v, u); + linop(m4).Tvmult_add(v, u); + + linop(m5).vmult(v, u); + linop(m5).Tvmult(v, u); + linop(m5).vmult_add(v, u); + linop(m5).Tvmult_add(v, u); + + linop(m6).vmult(v, u); + linop(m6).Tvmult(v, u); + linop(m6).vmult_add(v, u); + linop(m6).Tvmult_add(v, u); +} diff --git a/tests/lac/linear_operator_06.with_cxx11=on.output b/tests/lac/linear_operator_06.with_cxx11=on.output new file mode 100644 index 0000000000..b1e55fa541 --- /dev/null +++ b/tests/lac/linear_operator_06.with_cxx11=on.output @@ -0,0 +1,25 @@ + +DEAL::MyMatrix1::vmult +DEAL::MyMatrix1::Tvmult +DEAL::MyMatrix1::vmult +DEAL::MyMatrix1::Tvmult +DEAL::MyMatrix2::vmult +DEAL::MyMatrix2::Tvmult +DEAL::MyMatrix2::vmult_add +DEAL::MyMatrix2::Tvmult_add +DEAL::MyMatrix3::vmult +DEAL::MyMatrix3::Tvmult +DEAL::MyMatrix3::vmult_add +DEAL::MyMatrix3::Tvmult_add +DEAL::MyMatrix4::vmult +DEAL::MyMatrix4::Tvmult +DEAL::MyMatrix4::vmult_add +DEAL::MyMatrix4::Tvmult_add +DEAL::MyMatrix5::vmult +DEAL::MyMatrix5::Tvmult +DEAL::MyMatrix5::vmult_add +DEAL::MyMatrix5::Tvmult_add +DEAL::MyMatrix6::vmult +DEAL::MyMatrix6::Tvmult +DEAL::MyMatrix6::vmult_add +DEAL::MyMatrix6::Tvmult_add