]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Implement a LinearOperator class
authorMatthias Maier <tamiko@43-1.org>
Wed, 8 Apr 2015 10:13:36 +0000 (12:13 +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 [new file with mode: 0644]
tests/lac/linear_operator_01.cc [new file with mode: 0644]
tests/lac/linear_operator_01.with_cxx11=on.output [new file with mode: 0644]
tests/lac/linear_operator_02.cc [new file with mode: 0644]
tests/lac/linear_operator_02.with_cxx11=on.output [new file with mode: 0644]
tests/lac/linear_operator_03.cc [new file with mode: 0644]
tests/lac/linear_operator_03.with_cxx11=on.output [new file with mode: 0644]
tests/lac/linear_operator_04.cc [new file with mode: 0644]
tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output [new file with mode: 0644]

diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
new file mode 100644 (file)
index 0000000..e08f06c
--- /dev/null
@@ -0,0 +1,706 @@
+// ---------------------------------------------------------------------
+//
+// 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__linear_operator_h
+#define __deal2__linear_operator_h
+
+#include <deal.II/base/config.h>
+#include <deal.II/base/exceptions.h>
+#include <deal.II/lac/vector_memory.h>
+
+#include <functional>
+
+#ifdef DEAL_II_WITH_CXX11
+
+DEAL_II_NAMESPACE_OPEN
+
+// Forward declarations:
+
+template <typename Number> class Vector;
+
+template <typename Range, typename Domain> class LinearOperator;
+template <typename Range = Vector<double>,
+          typename Domain = Range,
+          typename Matrix>
+LinearOperator<Range, Domain> linop(const Matrix &matrix);
+
+/**
+ * A class to store the abstract concept of a linear operator.
+ *
+ * This is achieved by storing the concept of such operations in an
+ * "abstract" class LinearOperator that only holds knowledge of how to
+ * apply a linear operation via an std::function object vmult (vmult_add,
+ * Tvmult, Tvmult_add), as well as, information on how to create an
+ * appropriate Domain and Range vector.
+ *
+ * The primary purpose of this class is to provide syntactic sugar for
+ * complex matrix-vector operations and free the user from having to
+ * create, set up and handle intermediate storage locations by hand.
+ *
+ * As an example consider the operation $(A+k\,B)\,C$, where $A$, $B$ and
+ * $C$ denote (possible different) matrices. In order to construct a
+ * LinearOperator <code>op</code> that stores the knowledge of this
+ * operation, one can write:
+ *
+ * @code
+ * dealii::SparseMatrix<double> A, B, C;
+ * const double k = ...;
+ *
+ * // Setup and assembly of matrices
+ *
+ * const auto op_a = linop(A);
+ * const auto op_b = linop(B);
+ * const auto op_c = linop(C);
+ *
+ * const auto op = (op_a + k * op_b) * op_c;
+ * @endcode
+ *
+ * This class is only available if deal.II was configured with C++11
+ * support, i.e., if <code>DEAL_II_WITH_CXX11</code> is enabled during
+ * cmake configure.
+ *
+ * @author: Luca Heltai, Matthias Maier, 2015
+ */
+template <typename Range, typename Domain> class LinearOperator
+{
+public:
+
+  /**
+   * Create an empty LinearOperator object. All std::function objects are
+   * initialized with default variants that throw an exception upon
+   * invocation.
+   */
+  LinearOperator ()
+  {
+    vmult = [](Range &, const Domain &)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::vmult called"));
+    };
+
+    vmult_add = [](Range &, const Domain &)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::vmult_add called"));
+    };
+
+    Tvmult = [](Domain &, const Range &)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::Tvmult called"));
+    };
+
+    Tvmult_add = [](Domain &, const Range &)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::Tvmult_add called"));
+    };
+
+    reinit_range_vector = [](Range &, bool)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::reinit_range_vector method called"));
+    };
+
+    reinit_domain_vector = [](Domain &, bool)
+    {
+      Assert(false, ExcMessage("Uninitialized LinearOperator<Range, "
+                               "Domain>::reinit_domain_vector method called"));
+    };
+  }
+
+  /**
+   * Create a LinearOperator object from an object @p op for which the
+   * conversion function <code>linop</code> is defined.
+   */
+  template<typename Op>
+  LinearOperator (const Op &op)
+  {
+    *this = linop<Range, Domain, Op>(op);
+  }
+
+  /**
+   * Copy assignment operator for an object @p op for which the conversion
+   * function <code>linop</code> is defined.
+   */
+  template <typename Op>
+  LinearOperator<Range, Domain> &operator=(const Op &op)
+  {
+    return *this = linop<Range, Domain, Op>(op);
+  }
+
+  /**
+   * Application of the LinearOperator object to a vector @p u of the
+   * Domain space giving a vector @p v of the Range space.
+   */
+  std::function<const void(Range &v, const Domain &u)> vmult;
+
+  /**
+   * Application of the LinearOperator object to a vector @p u of the
+   * Domain space. The result is added to the vector @p v.
+   */
+  std::function<const void(Range &v, const Domain &u)> vmult_add;
+
+  /**
+   * Application of the transpose LinearOperator object to a vector @p u of
+   * the Range space giving a vector @p v of the Domain space.
+   */
+  std::function<const void(Domain &v, const Range &u)> Tvmult;
+
+  /**
+   * Application of the transpose LinearOperator object to a vector @p u of
+   * the Range space. The result is added to the vector @p v.
+   */
+  std::function<const void(Domain &v, const Range &u)> Tvmult_add;
+
+  /**
+   * Initializes a vector of the Range space to be directly usable as the
+   * destination parameter in an application of vmult. Similar to the
+   * reinit functions of the vector classes, the boolean determines whether
+   * a fast initalization is done, i.e., if it is set to false the content
+   * of the vector is set to 0.
+   */
+  std::function<void(Range &, bool)> reinit_range_vector;
+
+  /**
+   * Initializes a vector of the Domain space to be directly usable as the
+   * source parameter in an application of vmult. Similar to the reinit
+   * functions of the vector classes, the boolean determines whether a fast
+   * initalization is done, i.e., if it is set to false the content of the
+   * vector is set to 0.
+   */
+  std::function<void(Domain &, bool)> reinit_domain_vector;
+
+  /**
+   * A memory object for vectors of Range space used for intermediate
+   * storage during computations in vmult, etc.
+   */
+  mutable GrowingVectorMemory<Range> range_vector_memory;
+
+  /**
+   * A memory object for vectors of Domain space used for intermediate
+   * storage during computations in vmult, etc.
+   */
+  mutable GrowingVectorMemory<Domain> domain_vector_memory;
+
+
+  /**
+   * Addition with a LinearOperator @p second_op with the same Domain
+   * and Range.
+   */
+  LinearOperator<Range, Domain> &
+  operator+=(const LinearOperator<Range, Domain> &second_op)
+  {
+    *this = *this + second_op;
+    return *this;
+  }
+
+  /**
+   * Subtraction with a LinearOperator @p second_op with the same Domain
+   * and Range.
+   */
+  LinearOperator<Range, Domain> &
+  operator-=(const LinearOperator<Range, Domain> &second_op)
+  {
+    *this = *this - second_op;
+    return *this;
+  }
+
+  /**
+   * Concatenation of the LinearOperator with an endormorphism @p second_op
+   * on the Domain space.
+   */
+  LinearOperator<Range, Domain> &
+  operator*=(const LinearOperator<Domain, Domain> &second_op)
+  {
+    *this = *this * second_op;
+    return *this;
+  }
+};
+
+
+
+/**
+ * Addition of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}+\text{second_op})x:=\text{first_op}(x)+\text{second_op}(x)$
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator+(const LinearOperator<Range, Domain> &first_op,
+          const LinearOperator<Range, Domain> &second_op)
+{
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_range_vector = first_op.reinit_range_vector;
+  return_op.reinit_domain_vector = first_op.reinit_domain_vector;
+
+  // ensure to have valid computation objects by catching first_op and
+  // second_op by value
+
+  return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
+  {
+    first_op.vmult(v, u);
+    second_op.vmult_add(v, u);
+  };
+
+  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
+  {
+    first_op.vmult_add(v, u);
+    second_op.vmult_add(v, u);
+  };
+
+  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
+  {
+    second_op.Tvmult(v, u);
+    first_op.Tvmult_add(v, u);
+  };
+
+  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
+  {
+    second_op.Tvmult_add(v, u);
+    first_op.Tvmult_add(v, u);
+  };
+
+  return return_op;
+}
+
+
+/**
+ * Subtraction of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}-\text{second_op})x:=\text{first_op}(x)-\text{second_op}(x)$
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator-(const LinearOperator<Range, Domain> &first_op,
+          const LinearOperator<Range, Domain> &second_op)
+{
+  // implement with addition and scalar multiplication
+  return first_op + (-1. * second_op);
+}
+
+
+/**
+ * Concatenation of two linear operators @p first_op and @p second_op given
+ * by $(\text{first_op}*\text{second_op})x:=\text{first_op}(\text{second_op}(x))$
+ */
+template <typename Range, typename Intermediate, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Intermediate> &first_op,
+          const LinearOperator<Intermediate, Domain> &second_op)
+{
+  LinearOperator<Range, Domain> return_op;
+
+  return_op.reinit_domain_vector = second_op.reinit_domain_vector;
+  return_op.reinit_range_vector = first_op.reinit_range_vector;
+
+  // ensure to have valid computation objects by catching first_op and
+  // second_op by value
+
+  // Note: The range_vector_memory and domain_vector_memory objects inside
+  // the lambda exrepessions are the one belonging to the by-value captured
+  // objects first_op and second_op. Thus, they belong in essence to the
+  // std::function object vmult (vmult_add, ...) - and therefore have the
+  // same lifetime as the function object and not the LinearOperator
+  // parameters first_op and second_op.
+
+  return_op.vmult = [first_op, second_op](Range &v, const Domain &u)
+  {
+    Intermediate *i = second_op.range_vector_memory.alloc();
+    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+    second_op.vmult(*i, u);
+    first_op.vmult(v, *i);
+    second_op.range_vector_memory.free(i);
+  };
+
+  return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u)
+  {
+    Intermediate *i = second_op.range_vector_memory.alloc();
+    second_op.reinit_range_vector(*i, /*bool fast =*/ true);
+    second_op.vmult(*i, u);
+    first_op.vmult_add(v, *i);
+    second_op.range_vector_memory.free(i);
+  };
+
+  return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u)
+  {
+    Intermediate *i = first_op.domain_vector_memory.alloc();
+    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+    first_op.Tvmult(*i, u);
+    second_op.Tvmult(v, *i);
+    first_op.domain_vector_memory.free(i);
+  };
+
+  return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u)
+  {
+    Intermediate *i = first_op.domain_vector_memory.alloc();
+    first_op.reinit_domain_vector(*i, /*bool fast =*/ true);
+    first_op.Tvmult(*i, u);
+    second_op.Tvmult_add(v, *i);
+    first_op.domain_vector_memory.free(i);
+  };
+
+  return return_op;
+}
+
+
+/**
+ * Scalar multiplication of a ScalarOperator object from the left.
+ *
+ * The Domain and Range types must implement the following
+ * <code>operator*=</code> member functions accepting the appropriate
+ * scalar Number type for rescaling:
+ *
+ * @code
+ * Domain & operator *=(Number);
+ * Range & operator *=(Number);
+ * @endcode
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator*(typename Range::value_type  number,
+          const LinearOperator<Range, Domain> &op)
+{
+  static_assert(
+    std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
+    "Range and Domain must have implicitly convertible 'value_type's");
+
+  LinearOperator<Range, Domain> return_op = op;
+
+  // ensure to have valid computation objects by catching number and op by
+  // value
+
+  return_op.vmult = [number, op](Range &v, const Domain &u)
+  {
+    op.vmult(v,u);
+    v *= number;
+  };
+
+  return_op.vmult_add = [number, op](Range &v, const Domain &u)
+  {
+    v /= number;
+    op.vmult_add(v,u);
+    v *= number;
+  };
+
+  return_op.Tvmult = [number, op](Domain &v, const Range &u)
+  {
+    op.Tvmult(v,u);
+    v *= number;
+  };
+
+  return_op.Tvmult_add = [number, op](Domain &v, const Range &u)
+  {
+    v /= number;
+    op.Tvmult_add(v,u);
+    v *= number;
+  };
+
+  return return_op;
+}
+
+
+/**
+ * Scalar multiplication of a ScalarOperator object from the right.
+ *
+ * The Domain and Range types must implement the following
+ * <code>operator*=</code> member functions for rescaling:
+ *
+ * @code
+ * Domain & operator *=(Number);
+ * Range & operator *=(Number);
+ * @endcode
+ */
+template <typename Range, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Domain> &op,
+          typename Domain::value_type  number)
+{
+  static_assert(
+    std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
+    "Range and Domain must have implicitly convertible 'value_type's");
+
+  return number * op;
+}
+
+
+/**
+ * Returns the transpose linear operations of @p op.
+ */
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+transpose_linop(const LinearOperator<Range, Domain> &op)
+{
+  LinearOperator<Domain, Range> return_op;
+
+  return_op.reinit_range_vector = op.reinit_domain_vector;
+  return_op.reinit_domain_vector = op.reinit_range_vector;
+
+  return_op.vmult = op.Tvmult;
+  return_op.vmult_add = op.Tvmult_add;
+  return_op.Tvmult = op.vmult;
+  return_op.Tvmult_add = op.vmult_add;
+
+  return return_op;
+}
+
+
+/**
+ * Returns a LinearOperator that is the identity of the vector space
+ * Domain.
+ *
+ * The function takes an <code>std::function</code> object @p
+ * exemplar as an argument to initialize the reinit_range_vector
+ * and reinit_domain_vector objects of the LinearOperator object.
+ */
+template <typename Range>
+LinearOperator<Range, Range>
+identity_linop(const std::function<void(Range &, bool)> &exemplar)
+{
+  LinearOperator<Range, Range> return_op;
+
+  return_op.reinit_range_vector = exemplar;
+  return_op.reinit_domain_vector = exemplar;
+
+  return_op.vmult = [](Range &v, const Range &u)
+  {
+    v = u;
+  };
+
+  return_op.vmult_add = [](Range &v, const Range &u)
+  {
+    v += u;
+  };
+
+  return_op.Tvmult = [](Range &v, const Range &u)
+  {
+    v = u;
+  };
+
+  return_op.Tvmult_add = [](Range &v, const Range &u)
+  {
+    v += u;
+  };
+
+  return return_op;
+}
+
+
+/**
+ * Returns an object representing the inverse of the LinearOperator @p op.
+ *
+ * The function takes references @p solver and @p preconditioner to an
+ * iterative solver and a preconditioner that are used in the
+ * <code>vmult</code> and <code>Tvmult</code> implementations of the
+ * LinearOperator object.
+ *
+ * The LinearOperator object that is created stores a reference to @p
+ * solver and @p preconditioner. Thus, both objects must remain a valid
+ * reference for the whole lifetime of the LinearOperator object. Internal
+ * data structures of the @p solver object will be modified upon invocation
+ * of <code>vmult</code> or <code>Tvmult</code>.
+ */
+template <typename Solver, typename Preconditioner>
+LinearOperator<typename Solver::vector_type, typename Solver::vector_type>
+inverse_linop(Solver &solver,
+              const Preconditioner &preconditioner,
+              const LinearOperator<typename Solver::vector_type, typename Solver::vector_type> &op)
+{
+  typedef typename Solver::vector_type Vector;
+
+  LinearOperator<Vector, Vector> return_op;
+
+  return_op.reinit_range_vector = op.reinit_domain_vector;
+  return_op.reinit_domain_vector = op.reinit_range_vector;
+
+  return_op.vmult = [op, &solver, &preconditioner](Vector &v, const Vector &u)
+  {
+    solver.solve(op, v, u, preconditioner);
+  };
+
+  // Note: The range_vector_memory and domain_vector_memory objects inside
+  // the lambda exrepessions are the one belonging to the by-value captured
+  // objects first_op and second_op. Thus, they belong in essence to the
+  // std::function object vmult (vmult_add, ...) - and therefore have the
+  // same lifetime as the function object and not the LinearOperator
+  // parameter op.
+
+  return_op.vmult_add =
+    [op, &solver, &preconditioner](Vector &v, const Vector &u)
+  {
+    Vector *v2 = op.range_vector_memory.alloc();
+    op.reinit_range_vector(*v2, /*bool fast =*/ true);
+    solver.solve(op, *v2, u, preconditioner);
+    v += *v2;
+    op.range_vector_memory.free(v2);
+  };
+
+  return_op.Tvmult = [op, &solver, &preconditioner]( Vector &v, const
+                                                     Vector &u)
+  {
+    solver.solve(transpose_linop(op), v, u, preconditioner);
+  };
+
+  return_op.Tvmult =
+    [op, &solver, &preconditioner](Vector &v, const Vector &u)
+  {
+    Vector *v2 = op.range_vector_memory.alloc();
+    op.reinit_range_vector(*v2, /*bool fast =*/ true);
+    solver.solve(transpose_linop(op), *v2, u, preconditioner);
+    v += *v2;
+    op.range_vector_memory.free(v2);
+  };
+
+  return return_op;
+}
+
+
+/**
+ * A factory class that is responsible to create a reinit_range_vector
+ * object for a given pair of vector type Range and matrix type Matrix.
+ *
+ * The generic version of this class just calls Range::reinit() with the
+ * result of Matrix::m(). This class is specialized for more complicated
+ * data structures, such as TrilinosWrappers::MPI::Vector, etc.
+ */
+template<typename Range>
+class ReinitRangeFactory
+{
+public:
+
+  template <typename Matrix>
+  std::function<void(Range &, bool)>
+  operator()(const Matrix &matrix)
+  {
+    return [&matrix](Range &v, bool fast)
+    {
+      v.reinit(matrix.m(), fast);
+    };
+  }
+};
+
+
+/**
+ * A factory class that is responsible to create a reinit_domain_vector
+ * object for a given pair of vector type Domain and matrix type Matrix.
+ *
+ * The generic version of this class just calls Domain::reinit() with the
+ * result of Matrix::n(). This class is specialized for more complicated
+ * data structures, such as TrilinosWrappers::MPI::Vector, etc.
+ */
+template<typename Domain>
+class ReinitDomainFactory
+{
+public:
+
+  template <typename Matrix>
+  std::function<void(Domain &, bool)>
+  operator()(const Matrix &matrix)
+  {
+    return [&matrix](Domain &v, bool fast)
+    {
+      v.reinit(matrix.n(), fast);
+    };
+  }
+};
+
+
+/**
+ * A function that encapsulates generic @p matrix objects that act on a
+ * compatible Vector type into a LinearOperator. The LinearOperator object
+ * that is created stores a reference to the matrix object. Thus, @p matrix
+ * must remain a valid reference for the whole lifetime of the
+ * LinearOperator object.
+ *
+ * All changes made on @p matrix after the creation of the LinearOperator
+ * object are reflected by the operator object. For example, it is a valid
+ * procedure to first create a LinearOperator and resize, reassemble the
+ * matrix later.
+ *
+ * The Matrix class in question must provide the following minimal
+ * interface:
+ *
+ * @code
+ * class Matrix
+ * {
+ * public:
+ *   // (type specific) information how to create a Range and Domain vector
+ *   // with appropriate size and internal layout
+ *
+ *   // 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);
+ *
+ *   // 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.)
+ *   Tvmult_add(VECTOR &dst, const VECTOR &src);
+ * };
+ * @endcode
+ *
+ * @author: Matthias Maier, 2015
+ */
+template <typename Range = Vector<double>,
+          typename Domain = Range,
+          typename Matrix>
+LinearOperator<Range, Domain> linop(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...
+
+  return_op.reinit_range_vector =
+    ReinitRangeFactory<Range>().operator()(matrix);
+
+  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);
+  };
+
+  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);
+  };
+
+  return return_op;
+}
+
+
+DEAL_II_NAMESPACE_CLOSE
+
+#endif // DEAL_II_WITH_CXX11
+#endif
diff --git a/tests/lac/linear_operator_01.cc b/tests/lac/linear_operator_01.cc
new file mode 100644 (file)
index 0000000..4d601e3
--- /dev/null
@@ -0,0 +1,212 @@
+// ---------------------------------------------------------------------
+//
+// 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 the LinearOperator template on a trivial vector implementation
+// :: RightVector -> LeftVector
+
+#include "../tests.h"
+
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/vector_memory.templates.h>
+
+using namespace dealii;
+
+// Dummy vectors with different, non convertible types:
+
+struct LeftVector
+{
+  typedef double value_type;
+  double value;
+
+  LeftVector & operator *= (double scale)
+  {
+    value *= scale;
+    return *this;
+  }
+  LeftVector & operator /= (double scale)
+  {
+    value /= scale;
+    return *this;
+  }
+  LeftVector & operator += (const LeftVector &u)
+  {
+    value += u.value;
+    return *this;
+  }
+  int size() const
+  {
+    return 1;
+  }
+  std::size_t memory_consumption () const
+  {
+    return 1;
+  }
+};
+
+struct RightVector
+{
+  typedef double value_type;
+  double value;
+
+  RightVector & operator *= (double scale)
+  {
+    value *= scale;
+    return *this;
+  }
+  RightVector & operator /= (double scale)
+  {
+    value /= scale;
+    return *this;
+  }
+  RightVector & operator += (const RightVector &u)
+  {
+    value += u.value;
+    return *this;
+  }
+  int size() const
+  {
+    return 1;
+  }
+  std::size_t memory_consumption () const
+  {
+    return 1;
+  }
+};
+
+int main()
+{
+  initlog();
+
+  // Create to distinct linear operators:
+
+  LinearOperator<LeftVector, RightVector> multiply2;
+  multiply2.vmult = [](LeftVector &v, const RightVector &u)
+  {
+    v.value = 2 * u.value;
+  };
+  multiply2.vmult_add = [](LeftVector &v, const RightVector &u)
+  {
+    v.value += 2 * u.value;
+  };
+  multiply2.Tvmult = [](RightVector &v, const LeftVector &u)
+  {
+    v.value = 2 * u.value;
+  };
+  multiply2.Tvmult_add = [](RightVector &v, const LeftVector &u)
+  {
+    v.value += 2 * u.value;
+  };
+  multiply2.reinit_range_vector = [](LeftVector &, bool)
+  {
+  };
+  multiply2.reinit_domain_vector = [](RightVector &, bool)
+  {
+  };
+
+  auto multiply4 = multiply2;
+  multiply4.vmult = [](LeftVector &v, const RightVector &u)
+  {
+    v.value = 4 * u.value;
+  };
+  multiply4.vmult_add = [](LeftVector &v, const RightVector &u)
+  {
+    v.value += 4 * u.value;
+  };
+  multiply4.Tvmult = [](RightVector &v, const LeftVector &u)
+  {
+    v.value = 4 * u.value;
+  };
+  multiply4.Tvmult_add = [](RightVector &v, const LeftVector &u)
+  {
+    v.value += 4 * u.value;
+  };
+
+
+  // Small unit tests for all functions:
+
+  RightVector u = { 4. };
+  LeftVector v = { 0. };
+
+  // vmult, vmult_add
+
+  multiply2.vmult(v, u);
+  deallog << "2 * " << u.value << " = " << v.value << std::endl;
+
+  multiply4.vmult_add(v, u);
+  deallog << "... + 4 * " << u.value << " = " << v.value << std::endl;
+
+  multiply4.vmult(v, u);
+  deallog << "4 * " << u.value << " = " << v.value << std::endl;
+
+  multiply2.vmult_add(v, u);
+  deallog << "... + 2 * " << u.value << " = " << v.value << std::endl;
+
+  // Tvmult, Tvmult_add
+
+  v.value = 4.;
+
+  multiply2.Tvmult(u, v);
+  deallog << "2 * " << v.value << " = " << u.value << std::endl;
+
+  multiply4.Tvmult_add(u, v);
+  deallog << "... + 4 * " << v.value << " = " << u.value << std::endl;
+
+  multiply4.Tvmult(u, v);
+  deallog << "4 * " << v.value << " = " << u.value << std::endl;
+
+  multiply2.Tvmult_add(u, v);
+  deallog << "... + 2 * " << v.value << " = " << u.value << std::endl;
+
+  // operator+, operator-, operator+=, operator-=
+
+  auto test = multiply2 + multiply4;
+  test.vmult(v, u);
+  deallog << "(2 + 4) * " << u.value << " = " << v.value << std::endl;
+
+  test = multiply2 - multiply4;
+  test.vmult(v, u);
+  deallog << "(2 - 4) * " << u.value << " = " << v.value << std::endl;
+
+  test += multiply2;
+  test.vmult(v, u);
+  deallog << "(2 - 4 + 2) * " << u.value << " = " << v.value << std::endl;
+
+  test -= multiply4;
+  test.vmult(v, u);
+  deallog << "(2 - 4 + 2 - 4) * " << u.value << " = " << v.value << std::endl;
+
+  // operator* with scalar
+
+  test = 4. * multiply4;
+  test.vmult(v, u);
+  deallog << "(4 * 4) * " << u.value << " = " << v.value << std::endl;
+
+  test = multiply4 * 4.;
+  test.vmult(v, u);
+  deallog << "(4 * 4) * " << u.value << " = " << v.value << std::endl;
+
+  // operator* and transpose
+
+  auto test2 = transpose_linop(multiply2) * multiply4;
+  RightVector w = { 0. };
+  test2.vmult(w, u);
+  deallog << "(2 * 4) * " << u.value << " = " << w.value << std::endl;
+
+  // identity
+
+  auto test3 = identity_linop(test2.reinit_range_vector) + test2;
+  test3.vmult(w, u);
+  deallog << "(1 + 2 * 4) * " << u.value << " = " << w.value << std::endl;
+}
diff --git a/tests/lac/linear_operator_01.with_cxx11=on.output b/tests/lac/linear_operator_01.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..95b24c5
--- /dev/null
@@ -0,0 +1,17 @@
+
+DEAL::2 * 4.00000 = 8.00000
+DEAL::... + 4 * 4.00000 = 24.0000
+DEAL::4 * 4.00000 = 16.0000
+DEAL::... + 2 * 4.00000 = 24.0000
+DEAL::2 * 4.00000 = 8.00000
+DEAL::... + 4 * 4.00000 = 24.0000
+DEAL::4 * 4.00000 = 16.0000
+DEAL::... + 2 * 4.00000 = 24.0000
+DEAL::(2 + 4) * 24.0000 = 144.000
+DEAL::(2 - 4) * 24.0000 = -48.0000
+DEAL::(2 - 4 + 2) * 24.0000 = 0.00000
+DEAL::(2 - 4 + 2 - 4) * 24.0000 = -96.0000
+DEAL::(4 * 4) * 24.0000 = 384.000
+DEAL::(4 * 4) * 24.0000 = 384.000
+DEAL::(2 * 4) * 24.0000 = 192.000
+DEAL::(1 + 2 * 4) * 24.0000 = 216.000
diff --git a/tests/lac/linear_operator_02.cc b/tests/lac/linear_operator_02.cc
new file mode 100644 (file)
index 0000000..b2c34a1
--- /dev/null
@@ -0,0 +1,188 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Tests for the LinearOperator template with
+//   dealii::Vector<double>
+//   dealii::SparseMatrix<double>
+//   dealii::FullMatrix<double>
+
+#include "../tests.h"
+
+#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;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  static const int dim = 2;
+
+  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> a (sparsity_pattern);
+  SparseMatrix<double> b (sparsity_pattern);
+
+  QGauss<dim> quadrature(4);
+  MatrixCreator::create_laplace_matrix(mapping_q1, dof_handler, quadrature, a);
+  MatrixCreator::create_mass_matrix(mapping_q1, dof_handler, quadrature, b);
+
+
+  // Constructors and assignment:
+
+  auto op_a = linop(a);
+  auto op_b = linop(b);
+
+  {
+    LinearOperator<dealii::Vector<double>, dealii::Vector<double>> op_x (a);
+    op_a = a;
+    op_b = b;
+  }
+
+  // vmult:
+
+  Vector<double> u;
+  op_a.reinit_domain_vector(u, true);
+  for (unsigned int i = 0; i < u.size(); ++i) {
+    u[i] = (double)(i+1);
+  }
+
+  deallog << "u: " << u << std::endl;
+
+  Vector<double> v;
+  op_a.reinit_domain_vector(v, false);
+  Vector<double> w;
+  op_a.reinit_domain_vector(w, false);
+  Vector<double> x;
+  op_a.reinit_domain_vector(x, false);
+
+  op_a.vmult(v, u);
+  deallog << "Au: " << v << std::endl;
+
+  op_b.vmult(w, u);
+  deallog << "Bu: " << w << std::endl;
+
+  // operator+, operator-, operator+=, operator-=:
+
+  x = v;
+  x += w;
+  deallog << "Au+Bu: " << x << std::endl;
+
+  (op_a + op_b).vmult(x, u);
+  deallog << "(A+B)u: " << x << std::endl;
+
+  auto op_x = op_a;
+  op_x += op_b;
+  op_x.vmult(x, u);
+  deallog << "(A+=B)u: " << x << std::endl;
+
+  x = v;
+  x -= w;
+  deallog << "Au-Bu: " << x << std::endl;
+
+  (op_a - op_b).vmult(x, u);
+  deallog << "(A-B)u: " << x << std::endl;
+
+  op_x = op_a;
+  op_x -= op_b;
+  op_x.vmult(x, u);
+  deallog << "(A-=B)u: " << x << std::endl;
+
+  // operator*, operator*=
+
+  op_b.vmult(v,u);
+  op_a.vmult(w,v);
+  deallog << "(A(Bu)): " << w << std::endl;
+
+  (op_a * op_b).vmult(x, u);
+  deallog << "(A*B)u: " << x << std::endl;
+
+  op_x = op_a;
+  op_x *= op_b;
+  op_x.vmult(x, u);
+  deallog << "(A*=B)u: " << x << std::endl;
+
+  // solver interface:
+
+  SolverControl solver_control (1000, 1e-12);
+  SolverCG<> solver (solver_control);
+
+  deallog.depth_file(0);
+  solver.solve(op_b, v, u, PreconditionIdentity());
+  deallog.depth_file(3);
+  deallog << "solve(B, v, u): " << v << std::endl;
+
+  deallog.depth_file(0);
+  inverse_linop(solver, PreconditionIdentity(), op_b).vmult(v, u);
+  deallog.depth_file(3);
+  deallog << "inverse_linop(B)u: " << v << std::endl;
+
+  deallog.depth_file(0);
+  op_b.vmult(w, v);
+  deallog.depth_file(3);
+  deallog << "B(inverse_linop(B)u): " << w << std::endl;
+
+  deallog.depth_file(0);
+  (op_b * inverse_linop(solver, PreconditionIdentity(), op_b)).vmult(w, u);
+  deallog.depth_file(3);
+  deallog << "(B*inverse_linop(B))u: " << w << std::endl;
+
+  deallog.depth_file(0);
+  (inverse_linop(solver, PreconditionIdentity(), op_b) * op_b).vmult(w, u);
+  deallog.depth_file(3);
+  deallog << "(inverse_linop(B)*B)u: " << w << std::endl;
+
+  SolverControl inner_solver_control (1000, 1e-12);
+  SolverCG<> inner_solver (solver_control);
+
+  deallog.depth_file(0);
+  solver.solve(inverse_linop(inner_solver, PreconditionIdentity(), op_b), v, u,
+               PreconditionIdentity());
+  deallog.depth_file(3);
+  deallog << "solve(inverse_linop(B), v, u) == Bu: " << v << std::endl;
+
+  deallog.depth_file(0);
+  solver.solve(op_b + 0.005 * op_a, v, u, PreconditionIdentity());
+  deallog.depth_file(3);
+  deallog << "solve(B+0.5*0.01*A, v, u): " << v << std::endl;
+}
diff --git a/tests/lac/linear_operator_02.with_cxx11=on.output b/tests/lac/linear_operator_02.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..f9f1614
--- /dev/null
@@ -0,0 +1,39 @@
+
+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::Au: -1.500000000 -2.666666667 -2.000000000 -3.000000000 -2.333333333 -5.000000000 -3.500000000 -5.333333333 -9.333333333 0.5000000000 1.333333333 0.5000000000 1.166666667 -1.666666667 -1.666666667 2.000000000 6.000000000 3.000000000 1.000000000 3.000000000 1.666666667 9.000000000 4.000000000 3.333333333 1.500000000 
+DEAL::
+DEAL::Bu: 0.03125000000 0.09201388889 0.1145833333 0.2812500000 0.1788194444 0.4270833333 0.2552083333 0.5538194444 0.6631944444 0.3072916667 0.6753472222 0.1822916667 0.3923611111 0.8888888889 0.4878472222 0.4791666667 1.000000000 1.093750000 0.2864583333 0.5937500000 0.6371527778 1.281250000 0.6770833333 0.7170138889 0.3750000000 
+DEAL::
+DEAL::Au+Bu: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000 
+DEAL::
+DEAL::(A+B)u: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000 
+DEAL::
+DEAL::(A+=B)u: -1.468750000 -2.574652778 -1.885416667 -2.718750000 -2.154513889 -4.572916667 -3.244791667 -4.779513889 -8.670138889 0.8072916667 2.008680556 0.6822916667 1.559027778 -0.7777777778 -1.178819444 2.479166667 7.000000000 4.093750000 1.286458333 3.593750000 2.303819444 10.28125000 4.677083333 4.050347222 1.875000000 
+DEAL::
+DEAL::Au-Bu: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000 
+DEAL::
+DEAL::(A-B)u: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000 
+DEAL::
+DEAL::(A-=B)u: -1.531250000 -2.758680556 -2.114583333 -3.281250000 -2.512152778 -5.427083333 -3.755208333 -5.887152778 -9.996527778 0.1927083333 0.6579861111 0.3177083333 0.7743055556 -2.555555556 -2.154513889 1.520833333 5.000000000 1.906250000 0.7135416667 2.406250000 1.029513889 7.718750000 3.322916667 2.616319444 1.125000000 
+DEAL::
+DEAL::(A(Bu)): -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704 
+DEAL::
+DEAL::(A*B)u: -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704 
+DEAL::
+DEAL::(A*=B)u: -0.1073495370 -0.1866319444 -0.2039930556 -0.02199074074 -0.2893518519 -0.07465277778 -0.3703703704 0.03877314815 -0.2986111111 -0.1487268519 0.6250000000 -0.2201967593 -0.2123842593 0.4710648148 -0.4762731481 -0.1672453704 1.145833333 0.8049768519 -0.3211805556 -0.2199074074 -0.4939236111 1.570023148 -0.2034143519 -0.2300347222 -0.4094328704 
+DEAL::
+DEAL::solve(B, v, u): 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::inverse_linop(B)u: 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::B(inverse_linop(B)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*inverse_linop(B))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::(inverse_linop(B)*B)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::solve(inverse_linop(B), v, u) == Bu: 0.03125000000 0.09201388889 0.1145833333 0.2812500000 0.1788194444 0.4270833333 0.2552083333 0.5538194444 0.6631944444 0.3072916667 0.6753472222 0.1822916667 0.3923611111 0.8888888889 0.4878472222 0.4791666667 1.000000000 1.093750000 0.2864583333 0.5937500000 0.6371527778 1.281250000 0.6770833333 0.7170138889 0.3750000000 
+DEAL::
+DEAL::solve(B+0.5*0.01*A, v, u): 59.31059595 52.77843700 93.77361796 53.05371114 176.7862626 80.04635988 236.6925420 110.0293599 135.7027276 299.7370750 139.0997515 987.0719769 399.5852319 197.1836670 535.4087997 487.8703686 221.2510235 251.7282805 1564.703948 612.1971228 783.5848738 287.5180747 708.4784200 725.1386389 2062.398876 
+DEAL::
diff --git a/tests/lac/linear_operator_03.cc b/tests/lac/linear_operator_03.cc
new file mode 100644 (file)
index 0000000..85416d8
--- /dev/null
@@ -0,0 +1,208 @@
+// ---------------------------------------------------------------------
+//
+// 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.
+//
+// ---------------------------------------------------------------------
+
+// Tests for the LinearOperator template with
+//   dealii::BlockVector<double>
+//   dealii::BlockSparseMatrix<double>
+
+#include "../tests.h"
+
+#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/fe_system.h>
+#include <deal.II/fe/fe_values.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/block_sparse_matrix.h>
+#include <deal.II/lac/block_vector.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_gmres.h>
+#include <deal.II/lac/sparse_matrix.h>
+#include <deal.II/numerics/matrix_tools.h>
+
+#define PRINTME(name, var) \
+  deallog << name << ": [block 0] " << var.block(0) \
+                  << "  [block 1] " << var.block(1) << std::endl;
+
+
+using namespace dealii;
+
+int main()
+{
+  initlog();
+  deallog << std::setprecision(10);
+
+  static const int dim = 2;
+
+  Triangulation<dim> triangulation;
+  GridGenerator::hyper_cube (triangulation);
+  triangulation.refine_global(2);
+
+  MappingQ1<dim> mapping_q1;
+  FESystem<dim> fe(FE_Q<dim>(1), 1, FE_Q<dim>(1), 1);
+  DoFHandler<dim> dof_handler(triangulation);
+
+  dof_handler.distribute_dofs(fe);
+
+  std::vector<types::global_dof_index> dofs_per_component (2);
+  DoFTools::count_dofs_per_component (dof_handler, dofs_per_component);
+  const unsigned int n_u = dofs_per_component[0],
+                     n_p = dofs_per_component[1];
+
+  BlockDynamicSparsityPattern dsp(2, 2);
+  dsp.block(0, 0).reinit (n_u, n_u);
+  dsp.block(1, 0).reinit (n_p, n_u);
+  dsp.block(0, 1).reinit (n_u, n_p);
+  dsp.block(1, 1).reinit (n_p, n_p);
+  dsp.collect_sizes ();
+  DoFTools::make_sparsity_pattern (dof_handler, dsp);
+
+  BlockSparsityPattern sparsity_pattern;
+  sparsity_pattern.copy_from(dsp);
+  sparsity_pattern.compress();
+
+  BlockSparseMatrix<double> a (sparsity_pattern);
+  BlockSparseMatrix<double> b (sparsity_pattern);
+
+  for(unsigned int i = 0; i < a.n(); ++i)
+  {
+    a.set(i, i, 1.);
+    b.set(i, i, 5.);
+  }
+
+  // Constructors and assignment:
+
+  auto op_a = linop<BlockVector<double>>(a);
+  auto op_b = linop<BlockVector<double>>(b);
+
+  {
+    decltype(op_a) op_x (a);
+    op_a = a;
+    op_b = b;
+  }
+
+  // vmult:
+
+  BlockVector<double> u;
+  op_a.reinit_domain_vector(u, true);
+  for (unsigned int i = 0; i < u.size(); ++i) {
+    u[i] = (double)(i+1);
+  }
+
+  PRINTME("u", u);
+
+  BlockVector<double> v;
+  op_a.reinit_domain_vector(v, false);
+  BlockVector<double> w;
+  op_a.reinit_domain_vector(w, false);
+  BlockVector<double> x;
+  op_a.reinit_domain_vector(x, false);
+
+  op_a.vmult(v, u);
+  PRINTME("Au", v);
+
+  op_b.vmult(w, u);
+  PRINTME("Bu", w);
+
+  // operator+, operator-, operator+=, operator-=:
+
+  x = v;
+  x += w;
+  PRINTME("Au+Bu", x);
+
+  (op_a + op_b).vmult(x, u);
+  PRINTME("(A+B)u", x);
+
+  auto op_x = op_a;
+  op_x += op_b;
+  op_x.vmult(x, u);
+  PRINTME("(A+=B)u", x);
+
+  x = v;
+  x -= w;
+  PRINTME("Au-Bu", x);
+
+  (op_a - op_b).vmult(x, u);
+  PRINTME("(A-B)u", x);
+
+  op_x = op_a;
+  op_x -= op_b;
+  op_x.vmult(x, u);
+  PRINTME("(A-=B)u", x);
+
+  // operator*, operator*=
+
+  op_b.vmult(v,u);
+  op_a.vmult(w,v);
+  PRINTME("(A(Bu))", x);
+
+  (op_a * op_b).vmult(x, u);
+  PRINTME("(A*B)u", x);
+
+  op_x = op_a;
+  op_x *= op_b;
+  op_x.vmult(x, u);
+  PRINTME("(A*=B)u", x);
+
+  // solver interface:
+
+  SolverControl solver_control (1000, 1e-10);
+  SolverGMRES<BlockVector<double>> solver (solver_control);
+
+  deallog.depth_file(0);
+  solver.solve(op_b, v, u, PreconditionIdentity());
+  deallog.depth_file(3);
+  PRINTME("solve(B, v, u)", v);
+
+  deallog.depth_file(0);
+  inverse_linop(solver, PreconditionIdentity(), op_b).vmult(v, u);
+  deallog.depth_file(3);
+  PRINTME("inverse_linop(B)u", v);
+
+  deallog.depth_file(0);
+  op_b.vmult(w, v);
+  deallog.depth_file(3);
+  PRINTME("B(inverse_linop(B)u)", w);
+
+  deallog.depth_file(0);
+  (op_b * inverse_linop(solver, PreconditionIdentity(), op_b)).vmult(w, u);
+  deallog.depth_file(3);
+  PRINTME("(B*inverse_linop(B))u", w);
+
+  deallog.depth_file(0);
+  (inverse_linop(solver, PreconditionIdentity(), op_b) * op_b).vmult(w, u);
+  deallog.depth_file(3);
+  PRINTME("(inverse_linop(B)*B)u", w);
+
+  SolverControl inner_solver_control (1000, 1e-12);
+  SolverCG<BlockVector<double>> inner_solver (solver_control);
+
+  deallog.depth_file(0);
+  solver.solve(inverse_linop(inner_solver, PreconditionIdentity(), op_b), v, u,
+               PreconditionIdentity());
+  deallog.depth_file(3);
+  PRINTME("solve(inverse_linop(B), v, u) == Bu", v);
+
+  deallog.depth_file(0);
+  solver.solve(op_b + 0.005 * op_a, v, u, PreconditionIdentity());
+  deallog.depth_file(3);
+  PRINTME("solve(B+0.5*0.01*A, v, u)", v);
+}
diff --git a/tests/lac/linear_operator_03.with_cxx11=on.output b/tests/lac/linear_operator_03.with_cxx11=on.output
new file mode 100644 (file)
index 0000000..7d006da
--- /dev/null
@@ -0,0 +1,58 @@
+
+DEAL::u: [block 0] 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::  [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 
+DEAL::
+DEAL::Au: [block 0] 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::  [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 
+DEAL::
+DEAL::Bu: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000 
+DEAL::  [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000 
+DEAL::
+DEAL::Au+Bu: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000 
+DEAL::  [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000 
+DEAL::
+DEAL::(A+B)u: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000 
+DEAL::  [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000 
+DEAL::
+DEAL::(A+=B)u: [block 0] 6.000000000 12.00000000 18.00000000 24.00000000 30.00000000 36.00000000 42.00000000 48.00000000 54.00000000 60.00000000 66.00000000 72.00000000 78.00000000 84.00000000 90.00000000 96.00000000 102.0000000 108.0000000 114.0000000 120.0000000 126.0000000 132.0000000 138.0000000 144.0000000 150.0000000 
+DEAL::  [block 1] 156.0000000 162.0000000 168.0000000 174.0000000 180.0000000 186.0000000 192.0000000 198.0000000 204.0000000 210.0000000 216.0000000 222.0000000 228.0000000 234.0000000 240.0000000 246.0000000 252.0000000 258.0000000 264.0000000 270.0000000 276.0000000 282.0000000 288.0000000 294.0000000 300.0000000 
+DEAL::
+DEAL::Au-Bu: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000 
+DEAL::  [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000 
+DEAL::
+DEAL::(A-B)u: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000 
+DEAL::  [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000 
+DEAL::
+DEAL::(A-=B)u: [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000 
+DEAL::  [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000 
+DEAL::
+DEAL::(A(Bu)): [block 0] -4.000000000 -8.000000000 -12.00000000 -16.00000000 -20.00000000 -24.00000000 -28.00000000 -32.00000000 -36.00000000 -40.00000000 -44.00000000 -48.00000000 -52.00000000 -56.00000000 -60.00000000 -64.00000000 -68.00000000 -72.00000000 -76.00000000 -80.00000000 -84.00000000 -88.00000000 -92.00000000 -96.00000000 -100.0000000 
+DEAL::  [block 1] -104.0000000 -108.0000000 -112.0000000 -116.0000000 -120.0000000 -124.0000000 -128.0000000 -132.0000000 -136.0000000 -140.0000000 -144.0000000 -148.0000000 -152.0000000 -156.0000000 -160.0000000 -164.0000000 -168.0000000 -172.0000000 -176.0000000 -180.0000000 -184.0000000 -188.0000000 -192.0000000 -196.0000000 -200.0000000 
+DEAL::
+DEAL::(A*B)u: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000 
+DEAL::  [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000 
+DEAL::
+DEAL::(A*=B)u: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000 
+DEAL::  [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000 
+DEAL::
+DEAL::solve(B, v, u): [block 0] 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1.000000000 1.200000000 1.400000000 1.600000000 1.800000000 2.000000000 2.200000000 2.400000000 2.600000000 2.800000000 3.000000000 3.200000000 3.400000000 3.600000000 3.800000000 4.000000000 4.200000000 4.400000000 4.600000000 4.800000000 5.000000000 
+DEAL::  [block 1] 5.200000000 5.400000000 5.600000000 5.800000000 6.000000000 6.200000000 6.400000000 6.600000000 6.800000000 7.000000000 7.200000000 7.400000000 7.600000000 7.800000000 8.000000000 8.200000000 8.400000000 8.600000000 8.800000000 9.000000000 9.200000000 9.400000000 9.600000000 9.800000000 10.00000000 
+DEAL::
+DEAL::inverse_linop(B)u: [block 0] 0.2000000000 0.4000000000 0.6000000000 0.8000000000 1.000000000 1.200000000 1.400000000 1.600000000 1.800000000 2.000000000 2.200000000 2.400000000 2.600000000 2.800000000 3.000000000 3.200000000 3.400000000 3.600000000 3.800000000 4.000000000 4.200000000 4.400000000 4.600000000 4.800000000 5.000000000 
+DEAL::  [block 1] 5.200000000 5.400000000 5.600000000 5.800000000 6.000000000 6.200000000 6.400000000 6.600000000 6.800000000 7.000000000 7.200000000 7.400000000 7.600000000 7.800000000 8.000000000 8.200000000 8.400000000 8.600000000 8.800000000 9.000000000 9.200000000 9.400000000 9.600000000 9.800000000 10.00000000 
+DEAL::
+DEAL::B(inverse_linop(B)u): [block 0] 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::  [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 
+DEAL::
+DEAL::(B*inverse_linop(B))u: [block 0] 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::  [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 
+DEAL::
+DEAL::(inverse_linop(B)*B)u: [block 0] 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::  [block 1] 26.00000000 27.00000000 28.00000000 29.00000000 30.00000000 31.00000000 32.00000000 33.00000000 34.00000000 35.00000000 36.00000000 37.00000000 38.00000000 39.00000000 40.00000000 41.00000000 42.00000000 43.00000000 44.00000000 45.00000000 46.00000000 47.00000000 48.00000000 49.00000000 50.00000000 
+DEAL::
+DEAL::solve(inverse_linop(B), v, u) == Bu: [block 0] 5.000000000 10.00000000 15.00000000 20.00000000 25.00000000 30.00000000 35.00000000 40.00000000 45.00000000 50.00000000 55.00000000 60.00000000 65.00000000 70.00000000 75.00000000 80.00000000 85.00000000 90.00000000 95.00000000 100.0000000 105.0000000 110.0000000 115.0000000 120.0000000 125.0000000 
+DEAL::  [block 1] 130.0000000 135.0000000 140.0000000 145.0000000 150.0000000 155.0000000 160.0000000 165.0000000 170.0000000 175.0000000 180.0000000 185.0000000 190.0000000 195.0000000 200.0000000 205.0000000 210.0000000 215.0000000 220.0000000 225.0000000 230.0000000 235.0000000 240.0000000 245.0000000 250.0000000 
+DEAL::
+DEAL::solve(B+0.5*0.01*A, v, u): [block 0] 0.1998001998 0.3996003996 0.5994005994 0.7992007992 0.9990009990 1.198801199 1.398601399 1.598401598 1.798201798 1.998001998 2.197802198 2.397602398 2.597402597 2.797202797 2.997002997 3.196803197 3.396603397 3.596403596 3.796203796 3.996003996 4.195804196 4.395604396 4.595404595 4.795204795 4.995004995 
+DEAL::  [block 1] 5.194805195 5.394605395 5.594405594 5.794205794 5.994005994 6.193806194 6.393606394 6.593406593 6.793206793 6.993006993 7.192807193 7.392607393 7.592407592 7.792207792 7.992007992 8.191808192 8.391608392 8.591408591 8.791208791 8.991008991 9.190809191 9.390609391 9.590409590 9.790209790 9.990009990 
+DEAL::
diff --git a/tests/lac/linear_operator_04.cc b/tests/lac/linear_operator_04.cc
new file mode 100644 (file)
index 0000000..6f7f12d
--- /dev/null
@@ -0,0 +1,48 @@
+// ---------------------------------------------------------------------
+//
+// 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 that it is possible to instantiate a LinearOperator object for all
+// different kinds of Trilinos matrices and vectors
+
+#include "../tests.h"
+
+#include <deal.II/lac/linear_operator.h>
+#include <deal.II/lac/trilinos_block_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_vector.h>
+
+using namespace dealii;
+
+int main(int argc, char *argv[])
+{
+  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
+
+  initlog();
+  deallog << std::setprecision(10);
+
+  TrilinosWrappers::SparseMatrix a;
+
+  auto op_a  = linop<TrilinosWrappers::MPI::Vector>(a);
+  auto op_a2  = linop<TrilinosWrappers::Vector>(a);
+
+  TrilinosWrappers::BlockSparseMatrix b;
+
+  auto op_b = linop<TrilinosWrappers::MPI::BlockVector>(b);
+  auto op_b3 = linop<TrilinosWrappers::BlockVector>(b);
+
+  deallog << "OK" << std::endl;
+
+  return 0;
+}
diff --git a/tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output b/tests/lac/linear_operator_04.with_cxx11=on.with_trilinos=on.output
new file mode 100644 (file)
index 0000000..0fd8fc1
--- /dev/null
@@ -0,0 +1,2 @@
+
+DEAL::OK

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.