From f8eb97cf935cb81528f5aa63c35a1f768a7dc0d4 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Sun, 3 May 2015 15:49:03 +0200 Subject: [PATCH] rename and restructure --- include/deal.II/lac/computation.h | 760 -------------------------- include/deal.II/lac/linear_operator.h | 729 ++++++++++++++++++++++++ 2 files changed, 729 insertions(+), 760 deletions(-) delete mode 100644 include/deal.II/lac/computation.h diff --git a/include/deal.II/lac/computation.h b/include/deal.II/lac/computation.h deleted file mode 100644 index a583ef267d..0000000000 --- a/include/deal.II/lac/computation.h +++ /dev/null @@ -1,760 +0,0 @@ -// --------------------------------------------------------------------- -// -// Copyright (C) 2010 - 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. -// -// --------------------------------------------------------------------- - -#ifndef __deal2__computation_h -#define __deal2__computation_h - -#include -#include -#include - -#ifdef DEAL_II_WITH_CXX11 - -#include - -DEAL_II_NAMESPACE_OPEN - -// Forward declarations: - -template class LinearOperator; -template > class Computation; - - -/** - * A class to store a computation. - * - * The class essentially consists of std::function objects - * that store the knowledge of how to generate the result of a computation - * and store it in a vector: - * @code - * std::function result; - * std::function result_add; - * @endcode - * - * Similar to the LinearOperator class it also has knowledge about how to - * initialize a vector of the @p Range space: - * @code - * std::function reinit_vector; - * @endcode - * - * The primary purpose of this class is to allow lazy evaluation of - * expressions involving vectors and linear operators. As an example - * consider the addition of multiple vectors - * @code - * dealii::Vector a, b, c, d; - * // .. - * dealii::Vector result = a + b - c + d; - * @endcode - * or the computation of a residual $b-Ax$: - * @code - * dealii::SparseMatrix A; - * dealii::Vector b, x; - * // .. - * const auto op_a = linear_operator(A); - * - * const auto residual = b - op_a * x; - * @endcode - * The expression residual is of type - * Computation>. It stores references - * to A, b and x and defers the - * actual computation until result, or result_add - * are explicitly invoked, - * @code - * dealii::Vector y; - * residual.reinit_vector(y); - * residual.result(y); - * residual.result_add(y); - * @endcode - * or until the @p Computation object is implictly converted: - * @code - * dealii::Vector y; - * y = residual; - * y += residual; - * y -= residual; - * @endcode - * - * @author Matthias Maier, 2015 - * - * @ingroup LAOperators - */ -template class Computation -{ -public: - - /** - * Create an empty Computation object. All std::function - * member objects are initialized with default variants that throw an - * exception upon invocation. - */ - Computation () - { - result = [](Range &) - { - Assert(false, - ExcMessage("Uninitialized Computation::result called")); - }; - - result_add = [](Range &) - { - Assert(false, - ExcMessage("Uninitialized Computation::result_add called")); - }; - - reinit_vector = [](Range &, bool) - { - Assert(false, - ExcMessage("Uninitialized Computation::reinit_vector " - "method called")); - }; - } - - /** - * Default copy constructor. - */ - Computation (const Computation &) = default; - - /** - * Constructor that creates a Computation object from a reference vector - * @p u. The Computation returns @p u. - * - * The Computation object that is created stores a reference to @p u. - * Thus, the vector must remain a valid reference for the whole lifetime - * of the Computation object. All changes made on @p u after the creation - * of the Computation object are reflected by the operator object. - */ - Computation (const Range &u) - { - *this = u; - } - - /** - * Default copy assignment operator. - */ - Computation &operator=(const Computation &) = default; - - /** - * Copy assignment operator that creates a Computation object from a - * reference vector @p u. The Computation returns @p u. - * - * The Computation object that is created stores a reference to @p u. - * Thus, the vector must remain a valid reference for the whole lifetime - * of the Computation object. All changes made on @p u after the creation - * of the Computation object are reflected by the operator object. - */ - Computation &operator=(const Range &u) - { - result = [&u](Range &v) - { - v = u; - }; - - result_add = [&u](Range &v) - { - v += u; - }; - - reinit_vector = [&u](Range &v, bool fast) - { - v.reinit(u, fast); - }; - - return *this; - } - - /** - * Convert a Computation to its result. - * - * This conversion operator creates a vector of the Range space and calls - * result() on it. - */ - operator Range() const - { - Range result_vector; - - reinit_vector(result_vector, /*bool fast=*/ true); - result(result_vector); - - return result_vector; - } - - /** - * @name In-place vector space operations - */ - //@{ - - /** - * Addition with a Computation @p second_comp with the same @p Range. - */ - Computation &operator+=(const Computation &second_comp) - { - return *this = *this + second_comp; - } - - /** - * Subtraction with a Computation @p second_comp with the same @p Range. - */ - Computation &operator-=(const Computation &second_comp) - { - return *this = *this - second_comp; - } - - /** - * Add a constant @p offset (of the @p Range space) to the result of a - * computation. - */ - Computation &operator+=(const Range &offset) - { - return *this = *this + Computation(offset); - } - - /** - * Subract a constant @p offset (of the @p Range space) from the result - * of a computation. - */ - Computation &operator-=(const Range &offset) - { - return *this = *this - Computation(offset); - } - - /** - * Scalar multiplication of the Computation with a @p number. - */ - Computation &operator*=(typename Range::value_type number) - { - return *this = *this * number; - } - //@} - - /** - * Store the result of the computation in a vector v of the @p Range - * space. - */ - std::function result; - - /** - * Add the result of the computation to a vector v of the @p Range space. - */ - std::function result_add; - - /** - * Initializes a vector v of the Range space to be directly usable as the - * destination parameter in an application of result, or result_add. - * Similar to the reinit functions of the vector classes, the boolean - * determines whether a fast initialization is done, i.e., if it is set - * to false the content of the vector is set to 0. - */ - std::function reinit_vector; -}; - - -/** - * @name Vector space operations - */ -//@{ - -/** - * @relates Computation - * - * Addition of two Computation objects @p first_comp and @p second_comp given by - * vector space addition of the corresponding results. - * - * @ingroup LAOperators - */ -template -Computation -operator+(const Computation &first_comp, - const Computation &second_comp) -{ - Computation return_comp; - - return_comp.reinit_vector = first_comp.reinit_vector; - - // ensure to have valid computation objects by catching first_comp and - // second_comp by value - - return_comp.result = [first_comp, second_comp](Range &v) - { - first_comp.result(v); - second_comp.result_add(v); - }; - - return_comp.result_add = [first_comp, second_comp](Range &v) - { - first_comp.result_add(v); - second_comp.result_add(v); - }; - - return return_comp; -} - -/** - * @relates Computation - * - * Subtraction of two Computation objects @p first_comp and @p second_comp - * given by vector space addition of the corresponding results. - * - * @ingroup LAOperators - */ -template -Computation -operator-(const Computation &first_comp, - const Computation &second_comp) -{ - Computation return_comp; - - return_comp.reinit_vector = first_comp.reinit_vector; - - // ensure to have valid computation objects by catching first_comp and - // second_comp by value - - return_comp.result = [first_comp, second_comp](Range &v) - { - second_comp.result(v); - v *= -1.; - first_comp.result_add(v); - }; - - return_comp.result_add = [first_comp, second_comp](Range &v) - { - first_comp.result_add(v); - v *= -1.; - second_comp.result_add(v); - v *= -1.; - }; - - return return_comp; -} - -/** - * @relates Computation - * - * Scalar multiplication of a Computation objects @p comp with a scalar @p number - * given by a scaling Computation result with @p number. - * - * @ingroup LAOperators - */ -template -Computation -operator*(const Computation &comp, - typename Range::value_type number) -{ - Computation return_comp; - - return_comp.reinit_vector = comp.reinit_vector; - - return_comp.result = [comp, number](Range &v) - { - comp.result(v); - v *= number; - }; - - return_comp.result_add = [comp, number](Range &v) - { - v /= number; - comp.result_add(v); - v *= number; - }; - - return return_comp; -} - -/** - * @relates Computation - * - * Scalar multiplication of a Computation objects @p comp with a scalar @p number - * given by a scaling Computation result with @p number. - * - * @ingroup LAOperators - */ -template -Computation -operator*(typename Range::value_type number, - const Computation &comp) -{ - return comp * number; -} - -/** - * @relates Computation - * - * Add a constant @p offset (of the @p Range space) to the result of a - * Computation. - * - * @ingroup LAOperators - */ -template -Computation operator+(const Computation &comp, - const Range &offset) -{ - return comp + Computation(offset); -} - -/** - * @relates Computation - * - * Add a constant @p offset (of the @p Range space) to the result of a - * Computation. - * - * @ingroup LAOperators - */ -template -Computation operator+(const Range &offset, - const Computation &comp) -{ - return Computation(offset) + comp; -} - -/** - * @relates Computation - * - * Subtract a constant @p offset (of the @p Range space) from the result of a - * Computation. - * - * @ingroup LAOperators - */ -template -Computation operator-(const Computation &comp, - const Range &offset) -{ - return comp - Computation(offset); -} - - -/** - * @relates Computation - * - * Subtract a computational result from a constant @p offset (of the @p - * Range space). The result is a Computation object that applies this - * computation. - * - * @ingroup LAOperators - */ -template -Computation operator-(const Range &offset, - const Computation &comp) -{ - return Computation(offset) - comp; -} - -//@} - - -/** - * @name Creation of a Computation object - */ -//@{ - -namespace -{ - // Poor man's trait class that determines whether type T is a vector: - // FIXME: Implement this as a proprer type trait - similar to - // isBlockVector - - template - class has_vector_interface - { - template - static std::false_type test(...); - - template - static std::true_type test(decltype(&C::operator+=), - decltype(&C::operator-=), - decltype(&C::l2_norm)); - - 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, 0)) type; - }; -} - -/** - * @relates Computation - * - * Create a Computation object that stores the addition of two vectors. - * - * The Computation object that is created stores a reference to @p u and @p - * v. Thus, the vectors must remain valid references for the whole lifetime - * of the Computation object. All changes made on @p u or @p v after the - * creation of the Computation object are reflected by the operator object. - * - * @ingroup LAOperators - */ - -// FIXME: Only implement specialized variant for actual Vector types -template ::type::value>::type> -Computation operator+(const Range &u, const Range &v) -{ - Computation return_comp; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.reinit_vector = [&u](Range &x, bool fast) - { - x.reinit(u, fast); - }; - - return_comp.result = [&u, &v](Range &x) - { - x = u; - x += v; - }; - - return_comp.result_add = [&u, &v](Range &x) - { - x += u; - x += v; - }; - - return return_comp; -} - - -/** - * @relates Computation - * - * Create a Computation object that stores the addition of two vectors. - * - * The Computation object that is created stores a reference to @p u and @p - * v. Thus, the vectors must remain valid references for the whole lifetime - * of the Computation object. All changes made on @p u or @p v after the - * creation of the Computation object are reflected by the operator object. - * - * @ingroup LAOperators - */ - -// FIXME: Only implement specialized variant for actual Vector types -template ::type::value>::type> -Computation operator-(const Range &u, const Range &v) -{ - Computation return_comp; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.reinit_vector = [&u](Range &x, bool fast) - { - x.reinit(u, fast); - }; - - return_comp.result = [&u, &v](Range &x) - { - x = u; - x -= v; - }; - - return_comp.result_add = [&u, &v](Range &x) - { - x += u; - x -= v; - }; - - return return_comp; -} - - -/** - * @relates Computation - * - * Create a Computation object from a LinearOperator and a reference to a - * vector @p u of the Domain space. The object stores the computation - * $\text{op} \,u$ (in matrix notation). return - * (return_add) are implemented with vmult(__1,u) - * (vmult_add(__1,u)). - * - * The Computation object that is created stores a reference to @p u. - * Thus, the vector must remain a valid reference for the whole lifetime - * of the Computation object. All changes made on @p u after the creation - * of the Computation object are reflected by the operator object. - * - * @ingroup LAOperators - */ -template -Computation -operator*(const LinearOperator &op, - const Domain &u) -{ - Computation return_comp; - - return_comp.reinit_vector = op.reinit_range_vector; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.result = [op, &u](Range &v) - { - op.vmult(v, u); - }; - - return_comp.result_add = [op, &u](Range &v) - { - op.vmult_add(v, u); - }; - - return return_comp; -} - - -/** - * @relates Computation - * - * Create a Computation object from a LinearOperator and a reference to a - * vector @p u of the Range space. The object stores the computation - * $\text{op}^T \,u$ (in matrix notation). return - * (return_add) are implemented with Tvmult(__1,u) - * (Tvmult_add(__1,u)). - * - * The Computation object that is created stores a reference to @p u. - * Thus, the vector must remain a valid reference for the whole lifetime - * of the Computation object. All changes made on @p u after the creation - * of the Computation object are reflected by the operator object. - * - * @ingroup LAOperators - */ -template -Computation -operator*(const Range &u, - const LinearOperator &op) -{ - Computation return_comp; - - return_comp.reinit_vector = op.reinit_domain_vector; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.result = [op, &u](Domain &v) - { - op.Tvmult(v, u); - }; - - return_comp.result_add = [op, &u](Domain &v) - { - op.Tvmult_add(v, u); - }; - - return return_comp; -} - - -/** - * @relates Computation - * - * Composition of a Computation object with a LinearOperator. - * The object stores the computation $\text{op} \,comp$ (in matrix notation). - * - * @ingroup LAOperators - */ -template -Computation -operator*(const LinearOperator &op, - const Computation &comp) -{ - Computation return_comp; - - return_comp.reinit_vector = op.reinit_range_vector; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.result = [op, comp](Domain &v) - { - static GrowingVectorMemory vector_memory; - - Range *i = vector_memory.alloc(); - op.reinit_domain_vector(*i, /*bool fast =*/ true); - - comp.result(*i); - op.vmult(v, *i); - - vector_memory.free(i); - }; - - return_comp.result_add = [op, comp](Domain &v) - { - static GrowingVectorMemory vector_memory; - - Range *i = vector_memory.alloc(); - op.reinit_range_vector(*i, /*bool fast =*/ true); - - comp.result(*i); - op.vmult_add(v, *i); - - vector_memory.free(i); - }; - - return return_comp; -} - - -/** - * @relates Computation - * - * Composition of a Computation object with a LinearOperator. - * The object stores the computation $\text{op}^T \,comp$ (in matrix notation). - * - * @ingroup LAOperators - */ -template -Computation -operator*(const Computation &comp, - const LinearOperator &op) -{ - Computation return_comp; - - return_comp.reinit_vector = op.reinit_domain_vector; - - // ensure to have valid computation objects by catching op by value - // u is catched by reference - - return_comp.result = [op, comp](Domain &v) - { - static GrowingVectorMemory vector_memory; - - Range *i = vector_memory.alloc(); - op.reinit_range_vector(*i, /*bool fast =*/ true); - - comp.result(*i); - op.Tvmult(v, *i); - - vector_memory.free(i); - }; - - return_comp.result_add = [op, comp](Domain &v) - { - static GrowingVectorMemory vector_memory; - - Range *i = vector_memory.alloc(); - op.reinit_range_vector(*i, /*bool fast =*/ true); - - comp.result(*i); - op.Tvmult_add(v, *i); - - vector_memory.free(i); - }; - - return return_comp; -} - -//@} - -DEAL_II_NAMESPACE_CLOSE - -#endif // DEAL_II_WITH_CXX11 -#endif diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index 80ac97ed69..7eed984db9 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -47,6 +47,8 @@ template , typename Matrix> LinearOperator linear_operator (const Matrix &); +template > class PackagedOperation; + /** * A class to store the abstract concept of a linear operator. * @@ -1357,6 +1359,733 @@ block_diagonal_operator(const LinearOperatorstd::function objects + * that store the knowledge of how to generate the result of a computation + * and store it in a vector: + * @code + * std::function result; + * std::function result_add; + * @endcode + * + * Similar to the LinearOperator class it also has knowledge about how to + * initialize a vector of the @p Range space: + * @code + * std::function reinit_vector; + * @endcode + * + * The primary purpose of this class is to allow lazy evaluation of + * expressions involving vectors and linear operators. As an example + * consider the addition of multiple vectors + * @code + * dealii::Vector a, b, c, d; + * // .. + * dealii::Vector result = a + b - c + d; + * @endcode + * or the computation of a residual $b-Ax$: + * @code + * dealii::SparseMatrix A; + * dealii::Vector b, x; + * // .. + * const auto op_a = linear_operator(A); + * + * const auto residual = b - op_a * x; + * @endcode + * The expression residual is of type + * PackagedOperation>. It stores + * references to A, b and x and + * defers the actual computation until result, or + * result_add are explicitly invoked, + * @code + * dealii::Vector y; + * residual.reinit_vector(y); + * residual.result(y); + * residual.result_add(y); + * @endcode + * or until the @p PackagedOperation object is implictly converted: + * @code + * dealii::Vector y; + * y = residual; + * y += residual; + * y -= residual; + * @endcode + * + * @author Matthias Maier, 2015 + * + * @ingroup LAOperators + */ +template class PackagedOperation +{ +public: + + /** + * Create an empty PackagedOperation object. All std::function + * member objects are initialized with default variants that throw an + * exception upon invocation. + */ + PackagedOperation () + { + result = [](Range &) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::result called")); + }; + + result_add = [](Range &) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::result_add called")); + }; + + reinit_vector = [](Range &, bool) + { + Assert(false, + ExcMessage("Uninitialized PackagedOperation::reinit_vector " + "method called")); + }; + } + + /** + * Default copy constructor. + */ + PackagedOperation (const PackagedOperation &) = default; + + /** + * Constructor that creates a PackagedOperation object from a reference vector + * @p u. The PackagedOperation returns @p u. + * + * The PackagedOperation object that is created stores a reference to @p + * u. Thus, the vector must remain a valid reference for the whole + * lifetime of the PackagedOperation object. All changes made on @p u + * after the creation of the PackagedOperation object are reflected by + * the operator object. + */ + PackagedOperation (const Range &u) + { + *this = u; + } + + /** + * Default copy assignment operator. + */ + PackagedOperation &operator=(const PackagedOperation &) = default; + + /** Copy assignment operator that creates a PackagedOperation object from a + * reference vector @p u. The PackagedOperation returns @p u. + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime + * of the PackagedOperation object. All changes made on @p u after the creation + * of the PackagedOperation object are reflected by the operator object. + */ + PackagedOperation &operator=(const Range &u) + { + result = [&u](Range &v) + { + v = u; + }; + + result_add = [&u](Range &v) + { + v += u; + }; + + reinit_vector = [&u](Range &v, bool fast) + { + v.reinit(u, fast); + }; + + return *this; + } + + /** + * Convert a PackagedOperation to its result. + * + * This conversion operator creates a vector of the Range space and calls + * result() on it. + */ + operator Range() const + { + Range result_vector; + + reinit_vector(result_vector, /*bool fast=*/ true); + result(result_vector); + + return result_vector; + } + + /** + * @name In-place vector space operations + */ + //@{ + + /** + * Addition with a PackagedOperation @p second_comp with the same @p Range. + */ + PackagedOperation &operator+=(const PackagedOperation &second_comp) + { + return *this = *this + second_comp; + } + + /** + * Subtraction with a PackagedOperation @p second_comp with the same @p Range. + */ + PackagedOperation &operator-=(const PackagedOperation &second_comp) + { + return *this = *this - second_comp; + } + + /** + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + */ + PackagedOperation &operator+=(const Range &offset) + { + return *this = *this + PackagedOperation(offset); + } + + /** + * Subract a constant @p offset (of the @p Range space) from the result + * of a PackagedOperation. + */ + PackagedOperation &operator-=(const Range &offset) + { + return *this = *this - PackagedOperation(offset); + } + + /** + * Scalar multiplication of the PackagedOperation with a @p number. + */ + PackagedOperation &operator*=(typename Range::value_type number) + { + return *this = *this * number; + } + //@} + + /** + * Store the result of the PackagedOperation in a vector v of the @p Range + * space. + */ + std::function result; + + /** + * Add the result of the PackagedOperation to a vector v of the @p Range space. + */ + std::function result_add; + + /** + * Initializes a vector v of the Range space to be directly usable as the + * destination parameter in an application of result, or result_add. + * Similar to the reinit functions of the vector classes, the boolean + * determines whether a fast initialization is done, i.e., if it is set + * to false the content of the vector is set to 0. + */ + std::function reinit_vector; +}; + + +/** + * @name Vector space operations + */ +//@{ + +/** + * @relates PackagedOperation + * + * Addition of two PackagedOperation objects @p first_comp and + * @p second_comp given by vector space addition of the corresponding results. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator+(const PackagedOperation &first_comp, + const PackagedOperation &second_comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = first_comp.reinit_vector; + + // ensure to have valid PackagedOperation objects by catching first_comp and + // second_comp by value + + return_comp.result = [first_comp, second_comp](Range &v) + { + first_comp.result(v); + second_comp.result_add(v); + }; + + return_comp.result_add = [first_comp, second_comp](Range &v) + { + first_comp.result_add(v); + second_comp.result_add(v); + }; + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Subtraction of two PackagedOperation objects @p first_comp and @p second_comp + * given by vector space addition of the corresponding results. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator-(const PackagedOperation &first_comp, + const PackagedOperation &second_comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = first_comp.reinit_vector; + + // ensure to have valid PackagedOperation objects by catching first_comp and + // second_comp by value + + return_comp.result = [first_comp, second_comp](Range &v) + { + second_comp.result(v); + v *= -1.; + first_comp.result_add(v); + }; + + return_comp.result_add = [first_comp, second_comp](Range &v) + { + first_comp.result_add(v); + v *= -1.; + second_comp.result_add(v); + v *= -1.; + }; + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Scalar multiplication of a PackagedOperation objects @p comp with a + * scalar @p number given by a scaling PackagedOperation result with + * @p number. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const PackagedOperation &comp, + typename Range::value_type number) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = comp.reinit_vector; + + return_comp.result = [comp, number](Range &v) + { + comp.result(v); + v *= number; + }; + + return_comp.result_add = [comp, number](Range &v) + { + v /= number; + comp.result_add(v); + v *= number; + }; + + return return_comp; +} + +/** + * @relates PackagedOperation + * + * Scalar multiplication of a PackagedOperation objects @p comp with a + * scalar @p number given by a scaling PackagedOperation result with + * @p number. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(typename Range::value_type number, + const PackagedOperation &comp) +{ + return comp * number; +} + +/** + * @relates PackagedOperation + * + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator+(const PackagedOperation &comp, + const Range &offset) +{ + return comp + PackagedOperation(offset); +} + +/** + * @relates PackagedOperation + * + * Add a constant @p offset (of the @p Range space) to the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator+(const Range &offset, + const PackagedOperation &comp) +{ + return PackagedOperation(offset) + comp; +} + +/** + * @relates PackagedOperation + * + * Subtract a constant @p offset (of the @p Range space) from the result of a + * PackagedOperation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator-(const PackagedOperation &comp, + const Range &offset) +{ + return comp - PackagedOperation(offset); +} + + +/** + * @relates PackagedOperation + * + * Subtract a computational result from a constant @p offset (of the @p + * Range space). The result is a PackagedOperation object that applies this + * computation. + * + * @ingroup LAOperators + */ +template +PackagedOperation operator-(const Range &offset, + const PackagedOperation &comp) +{ + return PackagedOperation(offset) - comp; +} + +//@} + + +/** + * @name Creation of a PackagedOperation object + */ +//@{ + +namespace +{ + // Poor man's trait class that determines whether type T is a vector: + // FIXME: Implement this as a proprer type trait - similar to + // isBlockVector + + template + class has_vector_interface + { + template + static std::false_type test(...); + + template + static std::true_type test(decltype(&C::operator+=), + decltype(&C::operator-=), + decltype(&C::l2_norm)); + + 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, 0)) type; + }; +} + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the addition of two vectors. + * + * The PackagedOperation object that is created stores a reference to @p u and @p + * v. Thus, the vectors must remain valid references for the whole lifetime + * of the PackagedOperation object. All changes made on @p u or @p v after the + * creation of the PackagedOperation object are reflected by the operator object. + * + * @ingroup LAOperators + */ + +template ::type::value>::type> +PackagedOperation operator+(const Range &u, const Range &v) +{ + PackagedOperation return_comp; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.reinit_vector = [&u](Range &x, bool fast) + { + x.reinit(u, fast); + }; + + return_comp.result = [&u, &v](Range &x) + { + x = u; + x += v; + }; + + return_comp.result_add = [&u, &v](Range &x) + { + x += u; + x += v; + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object that stores the addition of two vectors. + * + * The PackagedOperation object that is created stores a reference to @p u + * and @p v. Thus, the vectors must remain valid references for the whole + * lifetime of the PackagedOperation object. All changes made on @p u or @p + * v after the creation of the PackagedOperation object are reflected by + * the operator object. + * + * @ingroup LAOperators + */ + +template ::type::value>::type> +PackagedOperation operator-(const Range &u, const Range &v) +{ + PackagedOperation return_comp; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.reinit_vector = [&u](Range &x, bool fast) + { + x.reinit(u, fast); + }; + + return_comp.result = [&u, &v](Range &x) + { + x = u; + x -= v; + }; + + return_comp.result_add = [&u, &v](Range &x) + { + x += u; + x -= v; + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object from a LinearOperator and a reference to a + * vector @p u of the Domain space. The object stores the PackagedOperation + * $\text{op} \,u$ (in matrix notation). return + * (return_add) are implemented with vmult(__1,u) + * (vmult_add(__1,u)). + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime of + * the PackagedOperation object. All changes made on @p u after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const LinearOperator &op, + const Domain &u) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_range_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.result = [op, &u](Range &v) + { + op.vmult(v, u); + }; + + return_comp.result_add = [op, &u](Range &v) + { + op.vmult_add(v, u); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Create a PackagedOperation object from a LinearOperator and a reference + * to a vector @p u of the Range space. The object stores the + * PackagedOperation $\text{op}^T \,u$ (in matrix notation). + * return (return_add) are implemented with + * Tvmult(__1,u) (Tvmult_add(__1,u)). + * + * The PackagedOperation object that is created stores a reference to @p u. + * Thus, the vector must remain a valid reference for the whole lifetime of + * the PackagedOperation object. All changes made on @p u after the + * creation of the PackagedOperation object are reflected by the operator + * object. + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const Range &u, + const LinearOperator &op) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_domain_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.result = [op, &u](Domain &v) + { + op.Tvmult(v, u); + }; + + return_comp.result_add = [op, &u](Domain &v) + { + op.Tvmult_add(v, u); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Composition of a PackagedOperation object with a LinearOperator. The + * object stores the computation $\text{op} \,comp$ (in matrix notation). + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const LinearOperator &op, + const PackagedOperation &comp) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_range_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.result = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_domain_vector(*i, /*bool fast =*/ true); + + comp.result(*i); + op.vmult(v, *i); + + vector_memory.free(i); + }; + + return_comp.result_add = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.result(*i); + op.vmult_add(v, *i); + + vector_memory.free(i); + }; + + return return_comp; +} + + +/** + * @relates PackagedOperation + * + * Composition of a PackagedOperation object with a LinearOperator. The + * object stores the computation $\text{op}^T \,comp$ (in matrix notation). + * + * @ingroup LAOperators + */ +template +PackagedOperation +operator*(const PackagedOperation &comp, + const LinearOperator &op) +{ + PackagedOperation return_comp; + + return_comp.reinit_vector = op.reinit_domain_vector; + + // ensure to have valid PackagedOperation objects by catching op by value + // u is catched by reference + + return_comp.result = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.result(*i); + op.Tvmult(v, *i); + + vector_memory.free(i); + }; + + return_comp.result_add = [op, comp](Domain &v) + { + static GrowingVectorMemory vector_memory; + + Range *i = vector_memory.alloc(); + op.reinit_range_vector(*i, /*bool fast =*/ true); + + comp.result(*i); + op.Tvmult_add(v, *i); + + vector_memory.free(i); + }; + + return return_comp; +} + +//@} + DEAL_II_NAMESPACE_CLOSE #endif // DEAL_II_WITH_CXX11 -- 2.39.5