]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Refactoring II: Reintroduce PackagedOperation and related tests
authorMatthias Maier <matthias.maier@iwr.uni-heidelberg.de>
Wed, 27 May 2015 12:23:37 +0000 (14:23 +0200)
committerMatthias Maier <tamiko@43-1.org>
Wed, 27 May 2015 20:08:10 +0000 (22:08 +0200)
This reintroduces PackagedOperation and related tests under a new header
file "packaged_operation.h". Tests renamed accordingly.

include/deal.II/lac/packaged_operation.h [new file with mode: 0644]
tests/lac/packaged_operation_01.cc [new file with mode: 0644]
tests/lac/packaged_operation_01.with_cxx11=on.output [new file with mode: 0644]
tests/lac/packaged_operation_02.cc [new file with mode: 0644]
tests/lac/packaged_operation_02.with_cxx11=on.output [new file with mode: 0644]

diff --git a/include/deal.II/lac/packaged_operation.h b/include/deal.II/lac/packaged_operation.h
new file mode 100644 (file)
index 0000000..13192cc
--- /dev/null
@@ -0,0 +1,830 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2014 - 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__packaged_operation_h
+#define __deal2__packaged_operation_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 Number> class Vector;
+template <typename Range, typename Domain> class LinearOperator;
+template <typename Range = Vector<double> > class PackagedOperation;
+
+
+/**
+ * A class to store a computation.
+ *
+ * The PackagedOperation class allows lazy evaluation of expressions
+ * involving vectors and linear operators. This is done by storing the
+ * computational expression and only performing the computation when either
+ * the object is implicitly converted to a vector object, or
+ * <code>apply</code> (or <code>apply_add</code>) is invoked by hand. This
+ * avoids unnecessary temporary storage of intermediate results.
+ *
+ * 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 &)> apply;
+ *   std::function<void(Range &)> apply_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
+ *
+ * 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);
+ *
+ *   dealii::Vector<double> 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>apply</code>, or
+ * <code>apply_add</code> are explicitly invoked,
+ * @code
+ *   dealii::Vector<double> y;
+ *   residual.reinit_vector(y);
+ *   residual.apply(y);
+ *   residual.apply_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 ()
+  {
+    apply = [](Range &)
+    {
+      Assert(false,
+             ExcMessage("Uninitialized PackagedOperation<Range>::apply called"));
+    };
+
+    apply_add = [](Range &)
+    {
+      Assert(false,
+             ExcMessage("Uninitialized PackagedOperation<Range>::apply_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)
+  {
+    apply = [&u](Range &v)
+    {
+      v = u;
+    };
+
+    apply_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>apply()</code> on it.
+   */
+  operator Range() const
+  {
+    Range result_vector;
+
+    reinit_vector(result_vector, /*bool fast=*/ true);
+    apply(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)> apply;
+
+  /**
+   * Add the result of the PackagedOperation to a vector v of the @p Range space.
+   */
+  std::function<void(Range &v)> apply_add;
+
+  /**
+   * Initializes a vector v of the Range space to be directly usable as the
+   * destination parameter in an application of apply, or apply_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.apply = [first_comp, second_comp](Range &v)
+  {
+    first_comp.apply(v);
+    second_comp.apply_add(v);
+  };
+
+  return_comp.apply_add = [first_comp, second_comp](Range &v)
+  {
+    first_comp.apply_add(v);
+    second_comp.apply_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.apply = [first_comp, second_comp](Range &v)
+  {
+    second_comp.apply(v);
+    v *= -1.;
+    first_comp.apply_add(v);
+  };
+
+  return_comp.apply_add = [first_comp, second_comp](Range &v)
+  {
+    first_comp.apply_add(v);
+    v *= -1.;
+    second_comp.apply_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;
+
+  // the trivial case: number is zero
+  if (number == 0.)
+    {
+      return_comp.apply = [](Range &v)
+      {
+        v = 0.;
+      };
+
+      return_comp.apply_add = [](Range &)
+      {
+      };
+    }
+  else
+    {
+      return_comp.apply = [comp, number](Range &v)
+      {
+        comp.apply(v);
+        v *= number;
+      };
+
+      return_comp.apply_add = [comp, number](Range &v)
+      {
+        v /= number;
+        comp.apply_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 proper 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.apply = [&u, &v](Range &x)
+  {
+    x = u;
+    x += v;
+  };
+
+  return_comp.apply_add = [&u, &v](Range &x)
+  {
+    x += u;
+    x += v;
+  };
+
+  return return_comp;
+}
+
+
+/**
+ * @relates PackagedOperation
+ *
+ * Create a PackagedOperation object that stores the subtraction 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.apply = [&u, &v](Range &x)
+  {
+    x = u;
+    x -= v;
+  };
+
+  return_comp.apply_add = [&u, &v](Range &x)
+  {
+    x += u;
+    x -= v;
+  };
+
+  return return_comp;
+}
+
+
+/**
+ * @relates PackagedOperation
+ *
+ * Create a PackagedOperation object that stores the scaling of a vector
+ * with a @p number.
+ *
+ * The PackagedOperation object that is created stores a reference to @p u.
+ * 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,
+                                   typename Range::value_type number)
+{
+  return PackagedOperation<Range>(u) * number;
+}
+
+
+/**
+ * @relates PackagedOperation
+ *
+ * Create a PackagedOperation object that stores the scaling of a vector
+ * with a @p number.
+ *
+ * The PackagedOperation object that is created stores a reference to @p u.
+ * 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*(typename Range::value_type number,
+                                   const Range &u)
+{
+  return number * PackagedOperation<Range>(u);
+}
+
+
+/**
+ * @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.apply = [op, &u](Range &v)
+  {
+    op.vmult(v, u);
+  };
+
+  return_comp.apply_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.apply = [op, &u](Domain &v)
+  {
+    op.Tvmult(v, u);
+  };
+
+  return_comp.apply_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.apply = [op, comp](Domain &v)
+  {
+    static GrowingVectorMemory<Range> vector_memory;
+
+    Range *i = vector_memory.alloc();
+    op.reinit_domain_vector(*i, /*bool fast =*/ true);
+
+    comp.apply(*i);
+    op.vmult(v, *i);
+
+    vector_memory.free(i);
+  };
+
+  return_comp.apply_add = [op, comp](Domain &v)
+  {
+    static GrowingVectorMemory<Range> vector_memory;
+
+    Range *i = vector_memory.alloc();
+    op.reinit_range_vector(*i, /*bool fast =*/ true);
+
+    comp.apply(*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.apply = [op, comp](Domain &v)
+  {
+    static GrowingVectorMemory<Range> vector_memory;
+
+    Range *i = vector_memory.alloc();
+    op.reinit_range_vector(*i, /*bool fast =*/ true);
+
+    comp.apply(*i);
+    op.Tvmult(v, *i);
+
+    vector_memory.free(i);
+  };
+
+  return_comp.apply_add = [op, comp](Domain &v)
+  {
+    static GrowingVectorMemory<Range> vector_memory;
+
+    Range *i = vector_memory.alloc();
+    op.reinit_range_vector(*i, /*bool fast =*/ true);
+
+    comp.apply(*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/tests/lac/packaged_operation_01.cc b/tests/lac/packaged_operation_01.cc
new file mode 100644 (file)
index 0000000..9acbc29
--- /dev/null
@@ -0,0 +1,84 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test PackagedOperation for
+//   dealii::Vector<double>
+//   dealii::SparseMatrix<double>
+
+#include "../tests.h"
+
+#include <deal.II/lac/vector.h>
+#include <deal.II/lac/packaged_operation.h>
+
+using namespace dealii;
+
+
+void test_applies(std::string description,
+                  const PackagedOperation<Vector<double> > &expr)
+{
+  // test apply
+  Vector<double> tmp = expr;
+  deallog << description << ": " << tmp << std::endl;
+
+  // test apply_add
+  for (auto &i : tmp)
+    i = 100.;
+  expr.apply_add(tmp);
+  deallog << "100. * 1_n + " << description << ": " << tmp << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  static const int dim = 2;
+
+  // Tests:
+
+  Vector<double> u(25);
+  for (unsigned int i = 0; i < u.size(); ++i) {
+    u[i] = (double)(i+1);
+  }
+
+  deallog << "u: " << u << std::endl;
+
+  // creation via operator+, operator-, operator*
+
+  test_applies("u + u", u + u);
+  test_applies("u - u", u - u);
+  test_applies("3. * u", 3. * u);
+  test_applies("u * 3.", u * 3.);
+
+  // creation via mixed operator+, operator-
+
+  auto expr = 2. * u;
+
+  test_applies("2. * u + u", expr + u);
+  test_applies("2. * u - u", expr - u);
+
+  test_applies("u + 2. * u", u + expr);
+  test_applies("u - 2. * u", u - expr);
+
+  // operator+, operator-, operator*
+
+  PackagedOperation<Vector<double> > expr2 = u;
+
+  test_applies("2. * u + u", expr + expr2);
+  test_applies("2. * u - u", expr - expr2);
+  test_applies("3. * u", 3. * expr2);
+  test_applies("u * 3.", expr2 * 3.);
+}
diff --git a/tests/lac/packaged_operation_01.with_cxx11=on.output b/tests/lac/packaged_operation_01.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..8efd013
--- /dev/null
@@ -0,0 +1,51 @@
+
+DEAL::u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::u + u: 2.000000000 4.000000000 6.000000000 8.000000000 10.00000000 12.00000000 14.00000000 16.00000000 18.00000000 20.00000000 22.00000000 24.00000000 26.00000000 28.00000000 30.00000000 32.00000000 34.00000000 36.00000000 38.00000000 40.00000000 42.00000000 44.00000000 46.00000000 48.00000000 50.00000000 
+DEAL::
+DEAL::100. * 1_n + u + u: 102.0000000 104.0000000 106.0000000 108.0000000 110.0000000 112.0000000 114.0000000 116.0000000 118.0000000 120.0000000 122.0000000 124.0000000 126.0000000 128.0000000 130.0000000 132.0000000 134.0000000 136.0000000 138.0000000 140.0000000 142.0000000 144.0000000 146.0000000 148.0000000 150.0000000 
+DEAL::
+DEAL::u - u: 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 0.000000000 
+DEAL::
+DEAL::100. * 1_n + u - u: 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 100.0000000 
+DEAL::
+DEAL::3. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + 3. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::u * 3.: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + u * 3.: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::2. * u + u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + 2. * u + u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::2. * u - u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + 2. * u - u: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::
+DEAL::u + 2. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + u + 2. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::u - 2. * u: -1.000000000 -2.000000000 -3.000000000 -4.000000000 -5.000000000 -6.000000000 -7.000000000 -8.000000000 -9.000000000 -10.00000000 -11.00000000 -12.00000000 -13.00000000 -14.00000000 -15.00000000 -16.00000000 -17.00000000 -18.00000000 -19.00000000 -20.00000000 -21.00000000 -22.00000000 -23.00000000 -24.00000000 -25.00000000 
+DEAL::
+DEAL::100. * 1_n + u - 2. * u: 99.00000000 98.00000000 97.00000000 96.00000000 95.00000000 94.00000000 93.00000000 92.00000000 91.00000000 90.00000000 89.00000000 88.00000000 87.00000000 86.00000000 85.00000000 84.00000000 83.00000000 82.00000000 81.00000000 80.00000000 79.00000000 78.00000000 77.00000000 76.00000000 75.00000000 
+DEAL::
+DEAL::2. * u + u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + 2. * u + u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::2. * u - u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + 2. * u - u: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::
+DEAL::3. * u: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + 3. * u: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
+DEAL::u * 3.: 3.000000000 6.000000000 9.000000000 12.00000000 15.00000000 18.00000000 21.00000000 24.00000000 27.00000000 30.00000000 33.00000000 36.00000000 39.00000000 42.00000000 45.00000000 48.00000000 51.00000000 54.00000000 57.00000000 60.00000000 63.00000000 66.00000000 69.00000000 72.00000000 75.00000000 
+DEAL::
+DEAL::100. * 1_n + u * 3.: 103.0000000 106.0000000 109.0000000 112.0000000 115.0000000 118.0000000 121.0000000 124.0000000 127.0000000 130.0000000 133.0000000 136.0000000 139.0000000 142.0000000 145.0000000 148.0000000 151.0000000 154.0000000 157.0000000 160.0000000 163.0000000 166.0000000 169.0000000 172.0000000 175.0000000 
+DEAL::
diff --git a/tests/lac/packaged_operation_02.cc b/tests/lac/packaged_operation_02.cc
new file mode 100644 (file)
index 0000000..070fa4e
--- /dev/null
@@ -0,0 +1,115 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2015 by the deal.II authors
+//
+// This file is part of the deal.II library.
+//
+// The deal.II library is free software; you can use it, redistribute
+// it, and/or modify it under the terms of the GNU Lesser General
+// Public License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+// The full text of the license can be found in the file LICENSE at
+// the top level of the deal.II distribution.
+//
+// ---------------------------------------------------------------------
+
+// Test PackagedOperation for
+//   dealii::SparseMatrix<double> with LinearOperator
+
+#include "../tests.h"
+
+#include <deal.II/lac/packaged_operation.h>
+
+// and a _lot_ of stuff to create a linera oprator
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/dofs/dof_handler.h>
+#include <deal.II/dofs/dof_tools.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/mapping_q.h>
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/tria.h>
+#include <deal.II/lac/dynamic_sparsity_pattern.h>
+#include <deal.II/lac/full_matrix.h>
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/precondition.h>
+#include <deal.II/lac/solver_cg.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+using namespace dealii;
+
+
+void test_applies(std::string description,
+                  const PackagedOperation<Vector<double> > &expr)
+{
+  // test apply
+  Vector<double> tmp = expr;
+  deallog << description << ": " << tmp << std::endl;
+
+  // test apply_add
+  for (auto &i : tmp)
+    i = 100.;
+  expr.apply_add(tmp);
+  deallog << "100. * 1_n + " << description << ": " << tmp << std::endl;
+}
+
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  static const int dim = 2;
+
+  // Create mass marix M, and an iterative inverse MInv:
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global(2);
+
+  MappingQ1<dim> mapping_q1;
+  FE_Q<dim> q1(1);
+  DoFHandler<dim> dof_handler(triangulation);
+  dof_handler.distribute_dofs(q1);
+
+  DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
+  DoFTools::make_sparsity_pattern(dof_handler, dsp);
+  dsp.compress();
+  SparsityPattern sparsity_pattern;
+  sparsity_pattern.copy_from(dsp);
+  sparsity_pattern.compress();
+
+  SparseMatrix<double> m (sparsity_pattern);
+
+  QGauss<dim> quadrature(4);
+  MatrixCreator::create_mass_matrix(mapping_q1, dof_handler, quadrature, m);
+
+  auto M = linear_operator(m);
+
+  SolverControl solver_control (1000, 1e-12);
+  SolverCG<> solver (solver_control);
+
+  auto MInv =  inverse_operator(M, solver, PreconditionIdentity());
+
+  // Tests:
+
+  Vector<double> u;
+  M.reinit_domain_vector(u, true);
+  for (unsigned int i = 0; i < u.size(); ++i) {
+    u[i] = (double)(i+1);
+  }
+
+  deallog << "u: " << u << std::endl;
+
+  deallog.depth_file(0);
+  Vector<double> b = MInv * u;
+  deallog.depth_file(3);
+  deallog << "b: " << b << std::endl;
+
+  test_applies("M * b", M * b);
+  test_applies("b * M", b * M);
+
+  auto expr = b;
+  test_applies("M * b", M * expr);
+  test_applies("b * M", expr * M);
+}
diff --git a/tests/lac/packaged_operation_02.with_cxx11=on.output b/tests/lac/packaged_operation_02.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..eb13437
--- /dev/null
@@ -0,0 +1,21 @@
+
+DEAL::u: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::b: 80.16326531 28.24489796 86.53061224 25.79591837 259.4285714 57.14285714 328.0000000 88.00000000 184.0000000 226.6122449 40.48979592 1446.693878 315.7551020 160.0000000 760.0000000 387.7551020 81.63265306 190.8571429 2286.693878 470.0408163 1140.571429 108.0816327 569.9591837 539.1020408 3018.448980 
+DEAL::
+DEAL::M * b: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + M * b: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::
+DEAL::b * M: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + b * M: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::
+DEAL::M * b: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + M * b: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::
+DEAL::b * M: 1.000000000 2.000000000 3.000000000 4.000000000 5.000000000 6.000000000 7.000000000 8.000000000 9.000000000 10.00000000 11.00000000 12.00000000 13.00000000 14.00000000 15.00000000 16.00000000 17.00000000 18.00000000 19.00000000 20.00000000 21.00000000 22.00000000 23.00000000 24.00000000 25.00000000 
+DEAL::
+DEAL::100. * 1_n + b * M: 101.0000000 102.0000000 103.0000000 104.0000000 105.0000000 106.0000000 107.0000000 108.0000000 109.0000000 110.0000000 111.0000000 112.0000000 113.0000000 114.0000000 115.0000000 116.0000000 117.0000000 118.0000000 119.0000000 120.0000000 121.0000000 122.0000000 123.0000000 124.0000000 125.0000000 
+DEAL::

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.