From: Matthias Maier <tamiko@43-1.org>
Date: Sun, 26 Apr 2015 21:36:41 +0000 (+0200)
Subject: [RFC] Add a mechanism to store a partially applied computation
X-Git-Tag: v8.3.0-rc1~179^2~15
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e06bc738b6cb1be050fd02f70ba4a579d9e72ae7;p=dealii.git

[RFC] Add a mechanism to store a partially applied computation
---

diff --git a/include/deal.II/lac/computation.h b/include/deal.II/lac/computation.h
new file mode 100644
index 0000000000..bcadf0fd75
--- /dev/null
+++ b/include/deal.II/lac/computation.h
@@ -0,0 +1,760 @@
+// ---------------------------------------------------------------------
+//
+// 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