+++ /dev/null
-// ---------------------------------------------------------------------
-//
-// 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 <deal.II/base/config.h>
-#include <deal.II/base/exceptions.h>
-#include <deal.II/lac/vector_memory.h>
-
-#ifdef DEAL_II_WITH_CXX11
-
-#include <functional>
-
-DEAL_II_NAMESPACE_OPEN
-
-// Forward declarations:
-
-template <typename Range, typename Domain> class LinearOperator;
-template <typename Range = Vector<double> > class Computation;
-
-
-/**
- * A class to store a computation.
- *
- * The class essentially consists of <code>std::function</code> objects
- * that store the knowledge of how to generate the result of a computation
- * and store it in a vector:
- * @code
- * std::function<void(Range &)> result;
- * std::function<void(Range &)> 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<void(Range &, bool)> 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<double> a, b, c, d;
- * // ..
- * dealii::Vector<double> result = a + b - c + d;
- * @endcode
- * or the computation of a residual $b-Ax$:
- * @code
- * dealii::SparseMatrix<double> A;
- * dealii::Vector<double> b, x;
- * // ..
- * const auto op_a = linear_operator(A);
- *
- * const auto residual = b - op_a * x;
- * @endcode
- * The expression <code>residual</code> is of type
- * <code>Computation<dealii::Vector<double>></code>. It stores references
- * to <code>A</code>, <code>b</code> and <code>x</code> and defers the
- * actual computation until <code>result</code>, or <code>result_add</code>
- * are explicitly invoked,
- * @code
- * dealii::Vector<double> 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<double> y;
- * y = residual;
- * y += residual;
- * y -= residual;
- * @endcode
- *
- * @author Matthias Maier, 2015
- *
- * @ingroup LAOperators
- */
-template <typename Range> class Computation
-{
-public:
-
- /**
- * Create an empty Computation object. All <code>std::function</code>
- * member objects are initialized with default variants that throw an
- * exception upon invocation.
- */
- Computation ()
- {
- result = [](Range &)
- {
- Assert(false,
- ExcMessage("Uninitialized Computation<Range>::result called"));
- };
-
- result_add = [](Range &)
- {
- Assert(false,
- ExcMessage("Uninitialized Computation<Range>::result_add called"));
- };
-
- reinit_vector = [](Range &, bool)
- {
- Assert(false,
- ExcMessage("Uninitialized Computation<Range>::reinit_vector "
- "method called"));
- };
- }
-
- /**
- * Default copy constructor.
- */
- Computation (const Computation<Range> &) = 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<Range> &operator=(const Computation<Range> &) = 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<Range> &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
- * <code>result()</code> 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<Range> &operator+=(const Computation<Range> &second_comp)
- {
- return *this = *this + second_comp;
- }
-
- /**
- * Subtraction with a Computation @p second_comp with the same @p Range.
- */
- Computation<Range> &operator-=(const Computation<Range> &second_comp)
- {
- return *this = *this - second_comp;
- }
-
- /**
- * Add a constant @p offset (of the @p Range space) to the result of a
- * computation.
- */
- Computation<Range> &operator+=(const Range &offset)
- {
- return *this = *this + Computation<Range>(offset);
- }
-
- /**
- * Subract a constant @p offset (of the @p Range space) from the result
- * of a computation.
- */
- Computation<Range> &operator-=(const Range &offset)
- {
- return *this = *this - Computation<Range>(offset);
- }
-
- /**
- * Scalar multiplication of the Computation with a @p number.
- */
- Computation<Range> &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<void(Range &v)> result;
-
- /**
- * Add the result of the computation to a vector v of the @p Range space.
- */
- std::function<void(Range &v)> 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<void(Range &v, bool fast)> 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 <typename Range>
-Computation<Range>
-operator+(const Computation<Range> &first_comp,
- const Computation<Range> &second_comp)
-{
- Computation<Range> 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 <typename Range>
-Computation<Range>
-operator-(const Computation<Range> &first_comp,
- const Computation<Range> &second_comp)
-{
- Computation<Range> 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 <typename Range>
-Computation<Range>
-operator*(const Computation<Range> &comp,
- typename Range::value_type number)
-{
- Computation<Range> 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 <typename Range>
-Computation<Range>
-operator*(typename Range::value_type number,
- const Computation<Range> &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 <typename Range>
-Computation<Range> operator+(const Computation<Range> &comp,
- const Range &offset)
-{
- return comp + Computation<Range>(offset);
-}
-
-/**
- * @relates Computation
- *
- * Add a constant @p offset (of the @p Range space) to the result of a
- * Computation.
- *
- * @ingroup LAOperators
- */
-template <typename Range>
-Computation<Range> operator+(const Range &offset,
- const Computation<Range> &comp)
-{
- return Computation<Range>(offset) + comp;
-}
-
-/**
- * @relates Computation
- *
- * Subtract a constant @p offset (of the @p Range space) from the result of a
- * Computation.
- *
- * @ingroup LAOperators
- */
-template <typename Range>
-Computation<Range> operator-(const Computation<Range> &comp,
- const Range &offset)
-{
- return comp - Computation<Range>(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 <typename Range>
-Computation<Range> operator-(const Range &offset,
- const Computation<Range> &comp)
-{
- return Computation<Range>(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 <typename T>
- class has_vector_interface
- {
- template <typename C>
- static std::false_type test(...);
-
- template <typename C>
- 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<T>(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 <typename Range,
- typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
-Computation<Range> operator+(const Range &u, const Range &v)
-{
- Computation<Range> 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 <typename Range,
- typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
-Computation<Range> operator-(const Range &u, const Range &v)
-{
- Computation<Range> 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). <code>return</code>
- * (<code>return_add</code>) are implemented with <code>vmult(__1,u)</code>
- * (<code>vmult_add(__1,u)</code>).
- *
- * 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 <typename Range, typename Domain>
-Computation<Range>
-operator*(const LinearOperator<Range, Domain> &op,
- const Domain &u)
-{
- Computation<Range> 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). <code>return</code>
- * (<code>return_add</code>) are implemented with <code>Tvmult(__1,u)</code>
- * (<code>Tvmult_add(__1,u)</code>).
- *
- * 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 <typename Range, typename Domain>
-Computation<Domain>
-operator*(const Range &u,
- const LinearOperator<Range, Domain> &op)
-{
- Computation<Range> 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 <typename Range, typename Domain>
-Computation<Range>
-operator*(const LinearOperator<Range, Domain> &op,
- const Computation<Domain> &comp)
-{
- Computation<Range> 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<Range> 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<Range> 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 <typename Range, typename Domain>
-Computation<Domain>
-operator*(const Computation<Range> &comp,
- const LinearOperator<Range, Domain> &op)
-{
- Computation<Range> 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<Range> 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<Range> 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
typename Matrix>
LinearOperator<Range, Domain> linear_operator (const Matrix &);
+template <typename Range = Vector<double> > class PackagedOperation;
+
/**
* A class to store the abstract concept of a linear operator.
*
//@}
+
+
+/**
+ * A class to store a computation.
+ *
+ * The class essentially consists of <code>std::function</code> objects
+ * that store the knowledge of how to generate the result of a computation
+ * and store it in a vector:
+ * @code
+ * std::function<void(Range &)> result;
+ * std::function<void(Range &)> 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<void(Range &, bool)> 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<double> a, b, c, d;
+ * // ..
+ * dealii::Vector<double> result = a + b - c + d;
+ * @endcode
+ * or the computation of a residual $b-Ax$:
+ * @code
+ * dealii::SparseMatrix<double> A;
+ * dealii::Vector<double> b, x;
+ * // ..
+ * const auto op_a = linear_operator(A);
+ *
+ * const auto residual = b - op_a * x;
+ * @endcode
+ * The expression <code>residual</code> is of type
+ * <code>PackagedOperation<dealii::Vector<double>></code>. It stores
+ * references to <code>A</code>, <code>b</code> and <code>x</code> and
+ * defers the actual computation until <code>result</code>, or
+ * <code>result_add</code> are explicitly invoked,
+ * @code
+ * dealii::Vector<double> 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<double> y;
+ * y = residual;
+ * y += residual;
+ * y -= residual;
+ * @endcode
+ *
+ * @author Matthias Maier, 2015
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range> class PackagedOperation
+{
+public:
+
+ /**
+ * Create an empty PackagedOperation object. All <code>std::function</code>
+ * member objects are initialized with default variants that throw an
+ * exception upon invocation.
+ */
+ PackagedOperation ()
+ {
+ result = [](Range &)
+ {
+ Assert(false,
+ ExcMessage("Uninitialized PackagedOperation<Range>::result called"));
+ };
+
+ result_add = [](Range &)
+ {
+ Assert(false,
+ ExcMessage("Uninitialized PackagedOperation<Range>::result_add called"));
+ };
+
+ reinit_vector = [](Range &, bool)
+ {
+ Assert(false,
+ ExcMessage("Uninitialized PackagedOperation<Range>::reinit_vector "
+ "method called"));
+ };
+ }
+
+ /**
+ * Default copy constructor.
+ */
+ PackagedOperation (const PackagedOperation<Range> &) = 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<Range> &operator=(const PackagedOperation<Range> &) = 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<Range> &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
+ * <code>result()</code> 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<Range> &operator+=(const PackagedOperation<Range> &second_comp)
+ {
+ return *this = *this + second_comp;
+ }
+
+ /**
+ * Subtraction with a PackagedOperation @p second_comp with the same @p Range.
+ */
+ PackagedOperation<Range> &operator-=(const PackagedOperation<Range> &second_comp)
+ {
+ return *this = *this - second_comp;
+ }
+
+ /**
+ * Add a constant @p offset (of the @p Range space) to the result of a
+ * PackagedOperation.
+ */
+ PackagedOperation<Range> &operator+=(const Range &offset)
+ {
+ return *this = *this + PackagedOperation<Range>(offset);
+ }
+
+ /**
+ * Subract a constant @p offset (of the @p Range space) from the result
+ * of a PackagedOperation.
+ */
+ PackagedOperation<Range> &operator-=(const Range &offset)
+ {
+ return *this = *this - PackagedOperation<Range>(offset);
+ }
+
+ /**
+ * Scalar multiplication of the PackagedOperation with a @p number.
+ */
+ PackagedOperation<Range> &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<void(Range &v)> result;
+
+ /**
+ * Add the result of the PackagedOperation to a vector v of the @p Range space.
+ */
+ std::function<void(Range &v)> 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<void(Range &v, bool fast)> 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 <typename Range>
+PackagedOperation<Range>
+operator+(const PackagedOperation<Range> &first_comp,
+ const PackagedOperation<Range> &second_comp)
+{
+ PackagedOperation<Range> 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 <typename Range>
+PackagedOperation<Range>
+operator-(const PackagedOperation<Range> &first_comp,
+ const PackagedOperation<Range> &second_comp)
+{
+ PackagedOperation<Range> 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 <typename Range>
+PackagedOperation<Range>
+operator*(const PackagedOperation<Range> &comp,
+ typename Range::value_type number)
+{
+ PackagedOperation<Range> 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 <typename Range>
+PackagedOperation<Range>
+operator*(typename Range::value_type number,
+ const PackagedOperation<Range> &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 <typename Range>
+PackagedOperation<Range> operator+(const PackagedOperation<Range> &comp,
+ const Range &offset)
+{
+ return comp + PackagedOperation<Range>(offset);
+}
+
+/**
+ * @relates PackagedOperation
+ *
+ * Add a constant @p offset (of the @p Range space) to the result of a
+ * PackagedOperation.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range>
+PackagedOperation<Range> operator+(const Range &offset,
+ const PackagedOperation<Range> &comp)
+{
+ return PackagedOperation<Range>(offset) + comp;
+}
+
+/**
+ * @relates PackagedOperation
+ *
+ * Subtract a constant @p offset (of the @p Range space) from the result of a
+ * PackagedOperation.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Range>
+PackagedOperation<Range> operator-(const PackagedOperation<Range> &comp,
+ const Range &offset)
+{
+ return comp - PackagedOperation<Range>(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 <typename Range>
+PackagedOperation<Range> operator-(const Range &offset,
+ const PackagedOperation<Range> &comp)
+{
+ return PackagedOperation<Range>(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 <typename T>
+ class has_vector_interface
+ {
+ template <typename C>
+ static std::false_type test(...);
+
+ template <typename C>
+ 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<T>(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 <typename Range,
+ typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
+PackagedOperation<Range> operator+(const Range &u, const Range &v)
+{
+ PackagedOperation<Range> 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 <typename Range,
+ typename = typename std::enable_if<has_vector_interface<Range>::type::value>::type>
+PackagedOperation<Range> operator-(const Range &u, const Range &v)
+{
+ PackagedOperation<Range> 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). <code>return</code>
+ * (<code>return_add</code>) are implemented with <code>vmult(__1,u)</code>
+ * (<code>vmult_add(__1,u)</code>).
+ *
+ * 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 <typename Range, typename Domain>
+PackagedOperation<Range>
+operator*(const LinearOperator<Range, Domain> &op,
+ const Domain &u)
+{
+ PackagedOperation<Range> 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).
+ * <code>return</code> (<code>return_add</code>) are implemented with
+ * <code>Tvmult(__1,u)</code> (<code>Tvmult_add(__1,u)</code>).
+ *
+ * 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 <typename Range, typename Domain>
+PackagedOperation<Domain>
+operator*(const Range &u,
+ const LinearOperator<Range, Domain> &op)
+{
+ PackagedOperation<Range> 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 <typename Range, typename Domain>
+PackagedOperation<Range>
+operator*(const LinearOperator<Range, Domain> &op,
+ const PackagedOperation<Domain> &comp)
+{
+ PackagedOperation<Range> 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<Range> 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<Range> 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 <typename Range, typename Domain>
+PackagedOperation<Domain>
+operator*(const PackagedOperation<Range> &comp,
+ const LinearOperator<Range, Domain> &op)
+{
+ PackagedOperation<Range> 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<Range> 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<Range> 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
projects / dealii.git / commitdiff