]> https://gitweb.dealii.org/ - dealii.git/commitdiff
rename and restructure
authorMatthias Maier <tamiko@43-1.org>
Sun, 3 May 2015 13:49:03 +0000 (15:49 +0200)
committerMatthias Maier <tamiko@43-1.org>
Mon, 11 May 2015 21:05:20 +0000 (23:05 +0200)
include/deal.II/lac/computation.h [deleted file]
include/deal.II/lac/linear_operator.h

diff --git a/include/deal.II/lac/computation.h b/include/deal.II/lac/computation.h
deleted file mode 100644 (file)
index a583ef2..0000000
+++ /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 <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
index 80ac97ed699851b5a829c757ff3f90b15570b9bb..7eed984db92627592eb95735f3ee160a3ecbb58e 100644 (file)
@@ -47,6 +47,8 @@ template <typename Range = Vector<double>,
           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.
  *
@@ -1357,6 +1359,733 @@ block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename
 
 //@}
 
+
+
+/**
+ * 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

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.