]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Generalize the "linop" constructor
authorMatthias Maier <tamiko@43-1.org>
Thu, 9 Apr 2015 00:40:54 +0000 (02:40 +0200)
committerMatthias Maier <tamiko@43-1.org>
Sun, 19 Apr 2015 20:55:06 +0000 (22:55 +0200)
include/deal.II/lac/linear_operator.h
tests/lac/linear_operator_06.cc [new file with mode: 0644]
tests/lac/linear_operator_06.with_cxx11=on.output [new file with mode: 0644]

index d95e21dcf2b2a05a701b8fabaaa1484a153e54b9..e437d304029c0c9da86c4e9251045bb7fa522773 100644 (file)
@@ -33,10 +33,17 @@ template <typename Number> class Vector;
 template <typename Number> class BlockVector;
 
 template <typename Range, typename Domain> class LinearOperator;
+
+template <typename Range = Vector<double>,
+          typename Domain = Range,
+          typename Exemplar,
+          typename Matrix>
+LinearOperator<Range, Domain> linop (const Exemplar &, const Matrix &);
+
 template <typename Range = Vector<double>,
           typename Domain = Range,
           typename Matrix>
-LinearOperator<Range, Domain> linop(const Matrix &matrix);
+LinearOperator<Range, Domain> linop (const Matrix &);
 
 /**
  * A class to store the abstract concept of a linear operator.
@@ -744,6 +751,120 @@ public:
 };
 
 
+namespace
+{
+  // A trait class that determines whether type T provides public
+  // (templated or non-templated) vmult_add and Tvmult_add member functions
+
+  template <typename Range, typename Domain, typename T>
+  class has_vmult_add
+  {
+    template <typename C>
+    static std::false_type test(...);
+
+    template <typename C>
+    static std::true_type test(decltype(&C::vmult_add),
+                               decltype(&C::Tvmult_add));
+
+    template <typename C>
+    static std::true_type test(decltype(&C::template vmult_add<Range>),
+                               decltype(&C::template Tvmult_add<Range>));
+
+    template <typename C>
+    static std::true_type test(decltype(&C::template vmult_add<Range, Domain>),
+                               decltype(&C::template Tvmult_add<Domain, Range>));
+
+  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)) type;
+  };
+
+
+  // A helper class to add the full matrix interface to a LinearOperator
+
+  template <typename Range, typename Domain>
+  class MatrixInterfaceWithVmultAdd
+  {
+  public:
+
+    template <typename Matrix>
+    void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+    {
+      op.vmult = [&matrix](Range &v, const Domain &u)
+      {
+        matrix.vmult(v,u);
+      };
+
+      op.vmult_add = [&matrix](Range &v, const Domain &u)
+      {
+        matrix.vmult_add(v,u);
+      };
+
+      op.Tvmult = [&matrix](Domain &v, const Range &u)
+      {
+        matrix.Tvmult(v,u);
+      };
+
+      op.Tvmult_add = [&matrix](Domain &v, const Range &u)
+      {
+        matrix.Tvmult_add(v,u);
+      };
+    }
+  };
+
+
+  // A helper class to add a reduced matrix interface to a LinearOperator
+  // (typically provided by Preconditioner classes)
+
+  template <typename Range, typename Domain>
+  class MatrixInterfaceWithoutVmultAdd
+  {
+  public:
+
+    template <typename Matrix>
+    void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+    {
+      op.vmult = [&matrix](Range &v, const Domain &u)
+      {
+        matrix.vmult(v,u);
+      };
+
+      op.vmult_add = [op](Range &v, const Domain &u)
+      {
+        // Capture op by value to have a valid GrowingVectorMemory object
+        // that is bound to the lifetime of this function object
+
+        Range *i = op.range_vector_memory.alloc();
+        op.reinit_range_vector(*i, /*bool fast =*/true);
+        op.vmult(*i, u);
+        v += *i;
+        op.range_vector_memory.free(i);
+      };
+
+      op.Tvmult = [&matrix](Domain &v, const Range &u)
+      {
+        matrix.Tvmult(v,u);
+      };
+
+      op.Tvmult_add = [op](Domain &v, const Range &u)
+      {
+        // Capture op by value to have a valid GrowingVectorMemory object
+        // that is bound to the lifetime of this function object
+
+        Domain *i = op.domain_vector_memory.alloc();
+        op.reinit_domain_vector(*i, /*bool fast =*/true);
+        op.Tvmult(*i, u);
+        v += *i;
+        op.domain_vector_memory.free(i);
+      };
+    }
+  };
+
+} /* namespace */
+
+
 /**
  * A function that encapsulates generic @p matrix objects that act on a
  * compatible Vector type into a LinearOperator. The LinearOperator object
@@ -769,60 +890,82 @@ public:
  *   // Application of matrix to vector src, writes the result into dst.
  *   vmult(VECTOR &dst, const VECTOR &src);
  *
- *   // Application of matrix to vector src, adds the result to dst.
- *   vmult_add(VECTOR &dst, const VECTOR &src);
- *
  *   // Application of the transpose of matrix to vector src, writes the
  *   // result into dst. (Depending on the usage of the linear operator
  *   // class this can be a dummy implementation throwing an error.)
  *   Tvmult(VECTOR &dst, const VECTOR &src);
+ * };
+ * @endcode
+ *
+ * The following (optional) interface is used if available:
+ *
+ * @code
+ * class Matrix
+ * {
+ * public:
+ *   // Application of matrix to vector src, adds the result to dst.
+ *   vmult_add(VECTOR &dst, const VECTOR &src);
  *
  *   // Application of the transpose of matrix to vector src, adds the
- *   // result to dst. (Depending on the usage of the linear operator class
- *   // this can be a dummy implementation throwing an error.)
+ *   // result to dst.
  *   Tvmult_add(VECTOR &dst, const VECTOR &src);
  * };
  * @endcode
  *
+ * If the Matrix does not provide <code>vmult_add</code> and
+ * <code>Tvmult_add</code>, they are implemented in terms of
+ * <code>vmult</code> and <code>Tvmult</code> (requiring intermediate
+ * storage).
+ *
  * @author: Matthias Maier, 2015
  */
 template <typename Range = Vector<double>,
           typename Domain = Range,
           typename Matrix>
 LinearOperator<Range, Domain> linop(const Matrix &matrix)
+{
+  // implement with the more generic variant below...
+  return linop<Range, Domain, Matrix, Matrix>(matrix, matrix);
+}
+
+
+/**
+ * Variant of above function that takes an Exemplar object @p exemplar as
+ * an additional reference. This object is used to populate the
+ * reinit_domain_vector and reinit_range_vector function objects. The
+ * reference @p matrix is used to construct vmult, Tvmult, etc.
+ *
+ * This variant can, for example, be used to encapsulate preconditioners
+ * (that typically do not expose any information about the underlying
+ * matrix).
+ *
+ * @author: Matthias Maier, 2015
+ */
+template <typename Range = Vector<double>,
+          typename Domain = Range,
+          typename Exemplar,
+          typename Matrix>
+LinearOperator<Range, Domain> linop(const Exemplar &exemplar,
+                                    const Matrix   &matrix)
 {
   LinearOperator<Range, Domain> return_op;
 
-  // Always store a reference to matrix in the lambda functions.
-  // This ensures that a modification of the matrix after the creation of a
-  // LinearOperator wrapper is respected - further a matrix cannot usually
-  // be copied...
+  // Always store a reference to matrix and exemplar in the lambda
+  // functions. This ensures that a modification of the matrix after the
+  // creation of a LinearOperator wrapper is respected - further a matrix
+  // or an exemplar cannot usually be copied...
 
   return_op.reinit_range_vector =
-    ReinitRangeFactory<Range>().operator()(matrix);
+    ReinitRangeFactory<Range>().operator()(exemplar);
 
   return_op.reinit_domain_vector =
-    ReinitDomainFactory<Domain>().operator()(matrix);
-
-  return_op.vmult = [&matrix](Range &v, const Domain &u)
-  {
-    matrix.vmult(v,u);
-  };
-
-  return_op.vmult_add = [&matrix](Range &v, const Domain &u)
-  {
-    matrix.vmult_add(v,u);
-  };
+    ReinitDomainFactory<Domain>().operator()(exemplar);
 
-  return_op.Tvmult = [&matrix](Domain &v, const Range &u)
-  {
-    matrix.Tvmult(v,u);
-  };
-
-  return_op.Tvmult_add = [&matrix](Domain &v, const Range &u)
-  {
-    matrix.Tvmult_add(v,u);
-  };
+  typename std::conditional<
+  has_vmult_add<Range, Domain, Matrix>::type::value,
+                MatrixInterfaceWithVmultAdd<Range, Domain>,
+                MatrixInterfaceWithoutVmultAdd<Range, Domain>>::type().
+                operator()(return_op, matrix);
 
   return return_op;
 }
diff --git a/tests/lac/linear_operator_06.cc b/tests/lac/linear_operator_06.cc
new file mode 100644 (file)
index 0000000..393b548
--- /dev/null
@@ -0,0 +1,239 @@
+// ---------------------------------------------------------------------
+//
+// 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 linop constructor for full and reduced interface:
+
+#include "../tests.h"
+
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/vector_memory.templates.h>
+#include <deal.II/lac/vector.h>
+
+class MyMatrix1
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  void vmult(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix1::vmult" << std::endl;
+  }
+
+  void Tvmult(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix1::Tvmult" << std::endl;
+  }
+
+};
+
+class MyMatrix2
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  void vmult(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix2::vmult" << std::endl;
+  }
+
+  void Tvmult(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix2::Tvmult" << std::endl;
+  }
+
+  void vmult_add(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix2::vmult_add" << std::endl;
+  }
+
+  void Tvmult_add(Vector<double> &, const Vector<double> &) const
+  {
+    deallog << "MyMatrix2::Tvmult_add" << std::endl;
+  }
+};
+
+class MyMatrix3
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  template<typename OutVector, typename InVector>
+  void vmult(OutVector &, const InVector &) const
+  {
+    deallog << "MyMatrix3::vmult" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void Tvmult(OutVector &, const InVector &) const
+  {
+    deallog << "MyMatrix3::Tvmult" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void vmult_add(OutVector &, const InVector &) const
+  {
+    deallog << "MyMatrix3::vmult_add" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void Tvmult_add(OutVector &, const InVector &) const
+  {
+    deallog << "MyMatrix3::Tvmult_add" << std::endl;
+  }
+};
+
+class MyMatrix4
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  template<typename OutVector, typename InVector>
+  void vmult(OutVector &, const InVector &, bool = true) const
+  {
+    deallog << "MyMatrix4::vmult" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void Tvmult(OutVector &, const InVector &, bool = true) const
+  {
+    deallog << "MyMatrix4::Tvmult" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void vmult_add(OutVector &, const InVector &, bool = true) const
+  {
+    deallog << "MyMatrix4::vmult_add" << std::endl;
+  }
+
+  template<typename OutVector, typename InVector>
+  void Tvmult_add(OutVector &, const InVector &, bool = true) const
+  {
+    deallog << "MyMatrix4::Tvmult_add" << std::endl;
+  }
+};
+
+class MyMatrix5
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  template<typename number>
+  void vmult(Vector<number> &, const Vector<number> &, bool = true) const
+  {
+    deallog << "MyMatrix5::vmult" << std::endl;
+  }
+
+  template<typename number>
+  void Tvmult(Vector<number> &, const Vector<number> &, bool = true) const
+  {
+    deallog << "MyMatrix5::Tvmult" << std::endl;
+  }
+
+  template<typename number>
+  void vmult_add(Vector<number> &, const Vector<number> &, bool = true) const
+  {
+    deallog << "MyMatrix5::vmult_add" << std::endl;
+  }
+
+  template<typename number>
+  void Tvmult_add(Vector<number> &, const Vector<number> &, bool = true) const
+  {
+    deallog << "MyMatrix5::Tvmult_add" << std::endl;
+  }
+};
+
+class MyMatrix6
+{
+public:
+  size_t m() const { return 1; };
+  size_t n() const { return 1; };
+
+  template<typename number, typename number2>
+  void vmult(Vector<number> &, const Vector<number2> &, bool = true) const
+  {
+    deallog << "MyMatrix6::vmult" << std::endl;
+  }
+
+  template<typename number, typename number2>
+  void Tvmult(Vector<number> &, const Vector<number2> &, bool = true) const
+  {
+    deallog << "MyMatrix6::Tvmult" << std::endl;
+  }
+
+  template<typename number, typename number2>
+  void vmult_add(Vector<number> &, const Vector<number2> &, bool = true) const
+  {
+    deallog << "MyMatrix6::vmult_add" << std::endl;
+  }
+
+  template<typename number, typename number2>
+  void Tvmult_add(Vector<number> &, const Vector<number2> &, bool = true) const
+  {
+    deallog << "MyMatrix6::Tvmult_add" << std::endl;
+  }
+};
+
+using namespace dealii;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  Vector<double> u (1);
+  Vector<double> v (1);
+
+  MyMatrix1 m1;
+  MyMatrix2 m2;
+  MyMatrix3 m3;
+  MyMatrix4 m4;
+  MyMatrix5 m5;
+  MyMatrix6 m6;
+
+  linop(m1).vmult(v, u);
+  linop(m1).Tvmult(v, u);
+  linop(m1).vmult_add(v, u);
+  linop(m1).Tvmult_add(v, u);
+
+  linop(m2).vmult(v, u);
+  linop(m2).Tvmult(v, u);
+  linop(m2).vmult_add(v, u);
+  linop(m2).Tvmult_add(v, u);
+
+  linop(m3).vmult(v, u);
+  linop(m3).Tvmult(v, u);
+  linop(m3).vmult_add(v, u);
+  linop(m3).Tvmult_add(v, u);
+
+  linop(m4).vmult(v, u);
+  linop(m4).Tvmult(v, u);
+  linop(m4).vmult_add(v, u);
+  linop(m4).Tvmult_add(v, u);
+
+  linop(m5).vmult(v, u);
+  linop(m5).Tvmult(v, u);
+  linop(m5).vmult_add(v, u);
+  linop(m5).Tvmult_add(v, u);
+
+  linop(m6).vmult(v, u);
+  linop(m6).Tvmult(v, u);
+  linop(m6).vmult_add(v, u);
+  linop(m6).Tvmult_add(v, u);
+}
diff --git a/tests/lac/linear_operator_06.with_cxx11=on.output b/tests/lac/linear_operator_06.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..b1e55fa
--- /dev/null
@@ -0,0 +1,25 @@
+
+DEAL::MyMatrix1::vmult
+DEAL::MyMatrix1::Tvmult
+DEAL::MyMatrix1::vmult
+DEAL::MyMatrix1::Tvmult
+DEAL::MyMatrix2::vmult
+DEAL::MyMatrix2::Tvmult
+DEAL::MyMatrix2::vmult_add
+DEAL::MyMatrix2::Tvmult_add
+DEAL::MyMatrix3::vmult
+DEAL::MyMatrix3::Tvmult
+DEAL::MyMatrix3::vmult_add
+DEAL::MyMatrix3::Tvmult_add
+DEAL::MyMatrix4::vmult
+DEAL::MyMatrix4::Tvmult
+DEAL::MyMatrix4::vmult_add
+DEAL::MyMatrix4::Tvmult_add
+DEAL::MyMatrix5::vmult
+DEAL::MyMatrix5::Tvmult
+DEAL::MyMatrix5::vmult_add
+DEAL::MyMatrix5::Tvmult_add
+DEAL::MyMatrix6::vmult
+DEAL::MyMatrix6::Tvmult
+DEAL::MyMatrix6::vmult_add
+DEAL::MyMatrix6::Tvmult_add

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.