From: Jean-Paul Pelteret <jppelteret@gmail.com>
Date: Mon, 16 Jan 2017 20:54:21 +0000 (+0100)
Subject: LinearOperators now derive from an arbitrary Payload class.
X-Git-Tag: v8.5.0-rc1~219^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=488628996e997ef7798369abc57c53c04d8a611c;p=dealii.git

LinearOperators now derive from an arbitrary Payload class.

This allows their functionality to be extended to support algebra class
other than deal.II's native classes.
---

diff --git a/doc/news/changes/minor/20170104Jean-PaulPelteretMatthiasMaier b/doc/news/changes/minor/20170104Jean-PaulPelteretMatthiasMaier
new file mode 100644
index 0000000000..af5c2c8a66
--- /dev/null
+++ b/doc/news/changes/minor/20170104Jean-PaulPelteretMatthiasMaier
@@ -0,0 +1,5 @@
+New: LinearOperator now derive from an arbitrary Payload class. This allows
+their functionality to be extended to support matrix and vector classes
+such as the wrappers to PETSc or Trilinos.
+<br>
+(Jean-Paul Pelteret, Matthias Maier, 2017/01/04)
diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h
index db8f6dcbf0..3ca3496061 100644
--- a/include/deal.II/lac/block_linear_operator.h
+++ b/include/deal.II/lac/block_linear_operator.h
@@ -27,35 +27,49 @@ DEAL_II_NAMESPACE_OPEN
 
 // Forward declarations:
 
+namespace internal
+{
+  namespace BlockLinearOperator
+  {
+    template<typename PayloadBlockType = internal::LinearOperator::EmptyPayload>
+    class EmptyBlockPayload;
+  }
+}
+
 template <typename Number> class BlockVector;
 
 template <typename Range = BlockVector<double>,
-          typename Domain = Range>
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
 class BlockLinearOperator;
 
 template <typename Range = BlockVector<double>,
           typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
           typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
 block_operator(const BlockMatrixType &matrix);
 
 template <size_t m, size_t n,
           typename Range = BlockVector<double>,
-          typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &);
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, n>, m> &);
 
 template <size_t m,
           typename Range = BlockVector<double>,
-          typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &);
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, m> &);
 
 template <size_t m,
           typename Range = BlockVector<double>,
-          typename Domain = Range>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op);
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> &op);
 
 // This is a workaround for a bug in <=gcc-4.7 that does not like partial
 // template default values in combination with local lambda expressions [1]
@@ -66,21 +80,24 @@ block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename
 
 template <typename Range = BlockVector<double>,
           typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>,
           typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
 block_diagonal_operator(const BlockMatrixType &block_matrix);
 
 template <typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Domain, Range>
-block_forward_substitution(const BlockLinearOperator<Range, Domain> &,
-                           const BlockLinearOperator<Domain, Range> &);
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_forward_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &,
+                           const BlockLinearOperator<Domain, Range, BlockPayload> &);
 
 template <typename Range = BlockVector<double>,
-          typename Domain = Range>
-LinearOperator<Domain, Range>
-block_back_substitution(const BlockLinearOperator<Range, Domain> &,
-                        const BlockLinearOperator<Domain, Range> &);
+          typename Domain = Range,
+          typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> >
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_back_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &,
+                        const BlockLinearOperator<Domain, Range, BlockPayload> &);
 
 // end of workaround
 
@@ -125,12 +142,12 @@ block_back_substitution(const BlockLinearOperator<Range, Domain> &,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-class BlockLinearOperator : public LinearOperator<Range, Domain>
+template <typename Range, typename Domain, typename BlockPayload>
+class BlockLinearOperator : public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
 {
 public:
 
-  typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType> BlockType;
+  typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> BlockType;
 
   /**
    * Create an empty BlockLinearOperator object.
@@ -139,10 +156,9 @@ public:
    * class LinearOperator are initialized with default variants that throw an
    * exception upon invocation.
    */
-  BlockLinearOperator()
-    : LinearOperator<Range, Domain>()
+  BlockLinearOperator(const BlockPayload &payload)
+    : LinearOperator<Range, Domain, typename BlockPayload::BlockType>(typename BlockPayload::BlockType(payload,payload))
   {
-
     n_block_rows = []() -> unsigned int
     {
       Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
@@ -165,7 +181,7 @@ public:
   /**
    * Default copy constructor.
    */
-  BlockLinearOperator(const BlockLinearOperator<Range, Domain> &) =
+  BlockLinearOperator(const BlockLinearOperator<Range, Domain, BlockPayload> &) =
     default;
 
   /**
@@ -176,7 +192,7 @@ public:
   template<typename Op>
   BlockLinearOperator(const Op &op)
   {
-    *this = block_operator<Range, Domain, Op>(op);
+    *this = block_operator<Range, Domain, BlockPayload, Op>(op);
   }
 
   /**
@@ -187,7 +203,7 @@ public:
   template<size_t m, size_t n>
   BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
   {
-    *this = block_operator<m, n, Range, Domain>(ops);
+    *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
   }
 
   /**
@@ -198,23 +214,23 @@ public:
   template<size_t m>
   BlockLinearOperator(const std::array<BlockType, m> &ops)
   {
-    *this = block_diagonal_operator<m, Range, Domain>(ops);
+    *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
   }
 
   /**
    * Default copy assignment operator.
    */
-  BlockLinearOperator<Range, Domain> &
-  operator=(const BlockLinearOperator<Range, Domain> &) = default;
+  BlockLinearOperator<Range, Domain, BlockPayload> &
+  operator=(const BlockLinearOperator<Range, Domain, BlockPayload> &) = default;
 
   /**
    * Templated copy assignment operator for an object @p op for which the
    * conversion function <code>block_operator</code> is defined.
    */
   template <typename Op>
-  BlockLinearOperator<Range, Domain> &operator=(const Op &op)
+  BlockLinearOperator<Range, Domain, BlockPayload> &operator=(const Op &op)
   {
-    *this = block_operator<Range, Domain, Op>(op);
+    *this = block_operator<Range, Domain, BlockPayload, Op>(op);
     return *this;
   }
 
@@ -224,10 +240,10 @@ public:
    * specialization.
    */
   template <size_t m, size_t n>
-  BlockLinearOperator<Range, Domain> &
+  BlockLinearOperator<Range, Domain, BlockPayload> &
   operator=(const std::array<std::array<BlockType, n>, m> &ops)
   {
-    *this = block_operator<m, n, Range, Domain>(ops);
+    *this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
     return *this;
   }
 
@@ -237,10 +253,10 @@ public:
    * operator calls the corresponding block_operator() specialization.
    */
   template <size_t m>
-  BlockLinearOperator<Range, Domain> &
+  BlockLinearOperator<Range, Domain, BlockPayload> &
   operator=(const std::array<BlockType, m> &ops)
   {
-    *this = block_diagonal_operator<m, Range, Domain>(ops);
+    *this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
     return *this;
   }
 
@@ -272,10 +288,10 @@ namespace internal
   {
     // Populate the LinearOperator interfaces with the help of the
     // BlockLinearOperator functions
-    template <typename Range, typename Domain>
+    template <typename Range, typename Domain, typename BlockPayload>
     inline void
     populate_linear_operator_functions(
-      dealii::BlockLinearOperator<Range, Domain> &op)
+      dealii::BlockLinearOperator<Range, Domain, BlockPayload> &op)
     {
       op.reinit_range_vector = [=](Range &v, bool omit_zeroing_entries)
       {
@@ -359,6 +375,42 @@ namespace internal
             op.block(j, i).Tvmult_add(v.block(i), u.block(j));
       };
     }
+
+
+    /**
+    * A dummy class for BlockLinearOperators that do not require any extensions
+    * to facilitate the operations of the block matrix or its subblocks.
+    *
+    * This is the Payload class typically associated with deal.II's native
+    * BlockSparseMatrix. To use Trilinos and PETSc BlockSparseMatrices it is
+    * necessary to initialize a BlockLinearOperator with their associated
+    * BlockPayload.
+    *
+    * @author Jean-Paul Pelteret, Matthias Maier, 2016
+    *
+    * @ingroup LAOperators
+    */
+    template<typename PayloadBlockType>
+    class EmptyBlockPayload
+    {
+    public:
+      /**
+      * Type of payload held by each subblock
+      */
+      typedef PayloadBlockType BlockType;
+
+      /**
+      * Default constructor
+      *
+      * Since this class does not do anything in particular and needs no special
+      * configuration, we have only one generic constructor that can be called
+      * under any conditions.
+      */
+      template <typename... Args>
+      EmptyBlockPayload (const Args &...)
+      { }
+    };
+
   } /*namespace BlockLinearOperator*/
 } /*namespace internal*/
 
@@ -382,13 +434,14 @@ namespace internal
  */
 template <typename Range,
           typename Domain,
+          typename BlockPayload,
           typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
 block_operator(const BlockMatrixType &block_matrix)
 {
-  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
 
-  BlockLinearOperator<Range, Domain> return_op;
+  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
 
   return_op.n_block_rows = [&block_matrix]() -> unsigned int
   {
@@ -445,16 +498,16 @@ block_operator(const BlockMatrixType &block_matrix)
  *
  * @ingroup LAOperators
  */
-template <size_t m, size_t n, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+template <size_t m, size_t n, typename Range, typename Domain, typename BlockPayload>
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, n>, m> &ops)
 {
   static_assert(m > 0 && n > 0,
                 "a blocked LinearOperator must consist of at least one block");
 
-  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
 
-  BlockLinearOperator<Range, Domain> return_op;
+  BlockLinearOperator<Range, Domain, BlockPayload> return_op ((BlockPayload())); // TODO: Create block payload so that this can be initialised correctly
 
   return_op.n_block_rows = []() -> unsigned int
   {
@@ -497,13 +550,14 @@ block_operator(const std::array<std::array<LinearOperator<typename Range::BlockT
  */
 template <typename Range,
           typename Domain,
+          typename BlockPayload,
           typename BlockMatrixType>
-BlockLinearOperator<Range, Domain>
+BlockLinearOperator<Range, Domain, BlockPayload>
 block_diagonal_operator(const BlockMatrixType &block_matrix)
 {
-  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
 
-  BlockLinearOperator<Range, Domain> return_op;
+  BlockLinearOperator<Range, Domain, BlockPayload> return_op (BlockPayload(block_matrix,block_matrix));
 
   return_op.n_block_rows = [&block_matrix]() -> unsigned int
   {
@@ -553,14 +607,14 @@ block_diagonal_operator(const BlockMatrixType &block_matrix)
  *
  * @ingroup LAOperators
  */
-template <size_t m, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
+template <size_t m, typename Range, typename Domain, typename BlockPayload>
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType>, m> &ops)
 {
   static_assert(m > 0,
                 "a blockdiagonal LinearOperator must consist of at least one block");
 
-  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
 
   std::array<std::array<BlockType, m>, m> new_ops;
 
@@ -597,15 +651,15 @@ block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockTyp
  *
  * @ingroup LAOperators
  */
-template <size_t m, typename Range, typename Domain>
-BlockLinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
+template <size_t m, typename Range, typename Domain, typename BlockPayload>
+BlockLinearOperator<Range, Domain, BlockPayload>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType> &op)
 {
   static_assert(m > 0,
                 "a blockdiagonal LinearOperator must consist of at least "
                 "one block");
 
-  typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+  typedef typename BlockLinearOperator<Range, Domain, BlockPayload>::BlockType BlockType;
   std::array<BlockType, m> new_ops;
   new_ops.fill(op);
 
@@ -655,12 +709,12 @@ block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-block_forward_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
-                           const BlockLinearOperator<Domain, Range> &diagonal_inverse)
+template <typename Range, typename Domain, typename BlockPayload>
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_forward_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
+                           const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
 {
-  LinearOperator<Range, Range> return_op;
+  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
 
   return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
   return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
@@ -767,12 +821,12 @@ block_forward_substitution(const BlockLinearOperator<Range, Domain> &block_opera
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-block_back_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
-                        const BlockLinearOperator<Domain, Range> &diagonal_inverse)
+template <typename Range, typename Domain, typename BlockPayload>
+LinearOperator<Domain, Range, typename BlockPayload::BlockType>
+block_back_substitution(const BlockLinearOperator<Range, Domain, BlockPayload> &block_operator,
+                        const BlockLinearOperator<Domain, Range, BlockPayload> &diagonal_inverse)
 {
-  LinearOperator<Range, Range> return_op;
+  LinearOperator<Range, Range, typename BlockPayload::BlockType> return_op ((typename BlockPayload::BlockType(diagonal_inverse)));
 
   return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector;
   return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector;
diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
index d563d93386..0fbfec3161 100644
--- a/include/deal.II/lac/linear_operator.h
+++ b/include/deal.II/lac/linear_operator.h
@@ -30,27 +30,40 @@ DEAL_II_NAMESPACE_OPEN
 
 // Forward declarations:
 
+namespace internal
+{
+  namespace LinearOperator
+  {
+    class EmptyPayload;
+  }
+}
+
 template <typename Number> class Vector;
 
-template <typename Range = Vector<double>, typename Domain = Range>
+template <typename Range = Vector<double>,
+          typename Domain = Range,
+          typename Payload = internal::LinearOperator::EmptyPayload>
 class LinearOperator;
 
 template <typename Range = Vector<double>,
           typename Domain = Range,
+          typename Payload = internal::LinearOperator::EmptyPayload,
           typename OperatorExemplar,
           typename Matrix>
-LinearOperator<Range, Domain> linear_operator (const OperatorExemplar &,
-                                               const Matrix &);
+LinearOperator<Range, Domain, Payload> linear_operator (const OperatorExemplar &,
+                                                        const Matrix &);
 
 template <typename Range = Vector<double>,
           typename Domain = Range,
+          typename Payload = internal::LinearOperator::EmptyPayload,
           typename Matrix>
-LinearOperator<Range, Domain> linear_operator (const Matrix &);
+LinearOperator<Range, Domain, Payload> linear_operator (const Matrix &);
 
 template <typename Range = Vector<double>,
-          typename Domain = Range>
-LinearOperator<Range, Domain>
-null_operator(const LinearOperator<Range, Domain> &);
+          typename Domain = Range,
+          typename Payload = internal::LinearOperator::EmptyPayload>
+LinearOperator<Range, Domain, Payload>
+null_operator(const LinearOperator<Range, Domain, Payload> &);
 
 
 /**
@@ -104,25 +117,58 @@ null_operator(const LinearOperator<Range, Domain> &);
  * this object to encapsulate matrix object of medium to large size (as a rule
  * of thumb, sparse matrices with a size $1000\times1000$, or larger).
  *
+ * @note In order to use Trilinos or PETSc sparse matrices and preconditioners
+ * in conjunction with the LinearOperator class, it is necessary to extend the
+ * functionality of the LinearOperator class by means of an additional Payload.
+ *
+ * For example, in order to construct an InverseOperator a call to a solver is
+ * required. Naturally these solvers don't have an interface to the
+ * LinearOperator (which, for example, may represent a composite operation).
+ * The TrilinosWrappers::internal::LinearOperator::TrilinosPayload therefore
+ * provides an interface extension to the LinearOperator so that it can be
+ * passed to the solver and used by the solver as if it were a Trilinos
+ * operator. This implies that all of the necessary functionality of the
+ * specific Trilinos operator has been overloaded within the Payload class.
+ * This includes operator-vector multiplication and inverse operator-vector
+ * multiplication, where the operator can be either a
+ * TrilinosWrappers::SparseMatrix or a TrilinosWrappers::PreconditionBase
+ * and the vector is a native Trilinos vector.
+ *
+ * Another case where payloads provide a crucial supplement to the
+ * LinearOperator class are when composite operations are constructed (via
+ * operator overloading). In this instance, it is again necessary to provide
+ * an interface that produces the result of this composite operation that is
+ * compatible with Trilinos operator used by Trilinos solvers.
+ *
+ * @note To ensure that the correct payload is provided, wrapper functions
+ * for linear operators have been provided within the respective
+ * TrilinosWrappers (and, in the future, PetscWrappers) namespaces.
+ *
  * @note 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
+ * @author Luca Heltai, Matthias Maier, 2015; Jean-Paul Pelteret, 2016
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain> class LinearOperator
+template <typename Range, typename Domain, typename Payload>
+class LinearOperator
+  : public Payload
 {
 public:
 
   /**
-   * Create an empty LinearOperator object. All <code>std::function</code>
-   * member objects are initialized with default variants that throw an
-   * exception upon invocation.
+   * Create an empty LinearOperator object.
+   * When a payload is passed to this constructor, the resulting operator is
+   * constructed with a functional payload.
+   * In either case, this constructor yields an object that can not actually
+   * be used for any linear operator operations, and will throw an exception
+   * upon invocation.
    */
-  LinearOperator()
+  LinearOperator(const Payload &payload = Payload())
     :
+    Payload (payload),
     is_null_operator(false)
   {
 
@@ -166,7 +212,7 @@ public:
   /**
    * Default copy constructor.
    */
-  LinearOperator (const LinearOperator<Range, Domain> &) = default;
+  LinearOperator (const LinearOperator<Range, Domain, Payload> &) = default;
 
   /**
    * Templated copy constructor that creates a LinearOperator object from an
@@ -174,27 +220,27 @@ public:
    * <code>linear_operator</code> is defined.
    */
   template<typename Op,
-           typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain>, Op>::value>::type>
+           typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>::type>
   LinearOperator (const Op &op)
   {
-    *this = linear_operator<Range, Domain, Op>(op);
+    *this = linear_operator<Range, Domain, Payload, Op>(op);
   }
 
   /**
    * Default copy assignment operator.
    */
-  LinearOperator<Range, Domain> &
-  operator=(const LinearOperator<Range, Domain> &) = default;
+  LinearOperator<Range, Domain, Payload> &
+  operator=(const LinearOperator<Range, Domain, Payload> &) = default;
 
   /**
    * Templated copy assignment operator for an object @p op for which the
    * conversion function <code>linear_operator</code> is defined.
    */
   template <typename Op,
-            typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain>, Op>::value>::type>
-  LinearOperator<Range, Domain> &operator=(const Op &op)
+            typename = typename std::enable_if<!std::is_base_of<LinearOperator<Range, Domain, Payload>, Op>::value>::type>
+  LinearOperator<Range, Domain, Payload> &operator=(const Op &op)
   {
-    *this = linear_operator<Range, Domain, Op>(op);
+    *this = linear_operator<Range, Domain, Payload, Op>(op);
     return *this;
   }
 
@@ -249,8 +295,8 @@ public:
    * Addition with a LinearOperator @p second_op with the same @p Domain and
    * @p Range.
    */
-  LinearOperator<Range, Domain> &
-  operator+=(const LinearOperator<Range, Domain> &second_op)
+  LinearOperator<Range, Domain, Payload> &
+  operator+=(const LinearOperator<Range, Domain, Payload> &second_op)
   {
     *this = *this + second_op;
     return *this;
@@ -260,8 +306,8 @@ public:
    * Subtraction with a LinearOperator @p second_op with the same @p Domain
    * and @p Range.
    */
-  LinearOperator<Range, Domain> &
-  operator-=(const LinearOperator<Range, Domain> &second_op)
+  LinearOperator<Range, Domain, Payload> &
+  operator-=(const LinearOperator<Range, Domain, Payload> &second_op)
   {
     *this = *this - second_op;
     return *this;
@@ -271,8 +317,8 @@ public:
    * Composition of the LinearOperator with an endomorphism @p second_op of
    * the @p Domain space.
    */
-  LinearOperator<Range, Domain> &
-  operator*=(const LinearOperator<Domain, Domain> &second_op)
+  LinearOperator<Range, Domain, Payload> &
+  operator*=(const LinearOperator<Domain, Domain, Payload> &second_op)
   {
     *this = *this * second_op;
     return *this;
@@ -282,7 +328,7 @@ public:
    * Scalar multiplication of the LinearOperator with @p number from the
    * right.
    */
-  LinearOperator<Range, Domain>
+  LinearOperator<Range, Domain, Payload>
   operator*=(typename Domain::value_type  number)
   {
     *this = *this * number;
@@ -314,10 +360,10 @@ public:
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator+(const LinearOperator<Range, Domain> &first_op,
-          const LinearOperator<Range, Domain> &second_op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator+(const LinearOperator<Range, Domain, Payload> &first_op,
+          const LinearOperator<Range, Domain, Payload> &second_op)
 {
   if (first_op.is_null_operator)
     {
@@ -329,7 +375,9 @@ operator+(const LinearOperator<Range, Domain> &first_op,
     }
   else
     {
-      LinearOperator<Range, Domain> return_op;
+      LinearOperator<Range, Domain, Payload> return_op (
+        static_cast<const Payload &>(first_op) + static_cast<const Payload &>(second_op)
+      );
 
       return_op.reinit_range_vector = first_op.reinit_range_vector;
       return_op.reinit_domain_vector = first_op.reinit_domain_vector;
@@ -375,14 +423,14 @@ operator+(const LinearOperator<Range, Domain> &first_op,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator-(const LinearOperator<Range, Domain> &first_op,
-          const LinearOperator<Range, Domain> &second_op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator-(const LinearOperator<Range, Domain, Payload> &first_op,
+          const LinearOperator<Range, Domain, Payload> &second_op)
 {
   if (first_op.is_null_operator)
     {
-      return -1 * second_op;
+      return -1. * second_op;
     }
   else if (second_op.is_null_operator)
     {
@@ -413,10 +461,10 @@ operator-(const LinearOperator<Range, Domain> &first_op,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
 operator*(typename Range::value_type  number,
-          const LinearOperator<Range, Domain> &op)
+          const LinearOperator<Range, Domain, Payload> &op)
 {
   static_assert(
     std::is_convertible<typename Range::value_type, typename Domain::value_type>::value,
@@ -432,7 +480,7 @@ operator*(typename Range::value_type  number,
     }
   else
     {
-      LinearOperator<Range, Domain> return_op = op;
+      LinearOperator<Range, Domain, Payload> return_op = op;
 
       // ensure to have valid computation objects by catching number and op by
       // value
@@ -483,9 +531,9 @@ operator*(typename Range::value_type  number,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Domain> &op,
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator*(const LinearOperator<Range, Domain, Payload> &op,
           typename Domain::value_type  number)
 {
   static_assert(
@@ -512,21 +560,23 @@ operator*(const LinearOperator<Range, Domain> &op,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Intermediate, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Intermediate> &first_op,
-          const LinearOperator<Intermediate, Domain> &second_op)
+template <typename Range, typename Intermediate, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+operator*(const LinearOperator<Range, Intermediate, Payload> &first_op,
+          const LinearOperator<Intermediate, Domain, Payload> &second_op)
 {
   if (first_op.is_null_operator || second_op.is_null_operator)
     {
-      LinearOperator<Range, Domain> return_op;
+      LinearOperator<Range, Domain, Payload> return_op;
       return_op.reinit_domain_vector = second_op.reinit_domain_vector;
       return_op.reinit_range_vector  = first_op.reinit_range_vector;
       return null_operator(return_op);
     }
   else
     {
-      LinearOperator<Range, Domain> return_op;
+      LinearOperator<Range, Domain, Payload> return_op (
+        static_cast<const Payload &>(first_op) * static_cast<const Payload &>(second_op)
+      );
 
       return_op.reinit_domain_vector = second_op.reinit_domain_vector;
       return_op.reinit_range_vector  = first_op.reinit_range_vector;
@@ -587,13 +637,17 @@ operator*(const LinearOperator<Range, Intermediate> &first_op,
  *
  * Return the transpose linear operations of @p op.
  *
+ * @author  Matthias Maier, 2015
+ *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-transpose_operator(const LinearOperator<Range, Domain> &op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Domain, Range, Payload>
+transpose_operator(const LinearOperator<Range, Domain, Payload> &op)
 {
-  LinearOperator<Domain, Range> return_op;
+  LinearOperator<Domain, Range, Payload> return_op (
+    op.transpose_payload()
+  );
 
   return_op.reinit_range_vector = op.reinit_domain_vector;
   return_op.reinit_domain_vector = op.reinit_range_vector;
@@ -606,6 +660,7 @@ transpose_operator(const LinearOperator<Range, Domain> &op)
   return return_op;
 }
 
+
 /**
  * @relates LinearOperator
  *
@@ -622,52 +677,55 @@ transpose_operator(const LinearOperator<Range, Domain> &op)
  * of the @p solver object will be modified upon invocation of
  * <code>vmult</code> or <code>Tvmult</code>.
  *
+ * @author Luca Heltai, Matthias Maier, Jean-Paul Pelteret, 2015
+ *
  * @ingroup LAOperators
  */
-template <typename Solver, typename Preconditioner>
-LinearOperator<typename Solver::vector_type, typename Solver::vector_type>
-inverse_operator(const LinearOperator<typename Solver::vector_type, typename Solver::vector_type> &op,
+template <typename Payload,
+          typename Solver, typename Preconditioner,
+          typename Range = typename Solver::vector_type, typename Domain = Range>
+LinearOperator<Domain, Range, Payload>
+inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
                  Solver &solver,
                  const Preconditioner &preconditioner)
 {
-  typedef typename Solver::vector_type Vector;
-
-  LinearOperator<Vector, Vector> return_op;
+  LinearOperator<Domain, Range, Payload> return_op (
+    op.inverse_payload(solver, preconditioner)
+  );
 
   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)
+  return_op.vmult = [op, &solver, &preconditioner](Range &v, const Domain &u)
   {
     op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/ false);
     solver.solve(op, v, u, preconditioner);
   };
 
   return_op.vmult_add =
-    [op, &solver, &preconditioner](Vector &v, const Vector &u)
+    [op, &solver, &preconditioner](Range &v, const Domain &u)
   {
-    static GrowingVectorMemory<typename Solver::vector_type> vector_memory;
+    static GrowingVectorMemory<Range> vector_memory;
 
-    Vector *v2 = vector_memory.alloc();
+    Range *v2 = vector_memory.alloc();
     op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/ false);
     solver.solve(op, *v2, u, preconditioner);
     v += *v2;
     vector_memory.free(v2);
   };
 
-  return_op.Tvmult = [op, &solver, &preconditioner]( Vector &v, const
-                                                     Vector &u)
+  return_op.Tvmult = [op, &solver, &preconditioner](Range &v, const Domain &u)
   {
     op.reinit_range_vector(v, /*bool omit_zeroing_entries =*/ false);
     solver.solve(transpose_operator(op), v, u, preconditioner);
   };
 
   return_op.Tvmult_add =
-    [op, &solver, &preconditioner](Vector &v, const Vector &u)
+    [op, &solver, &preconditioner](Range &v, const Domain &u)
   {
-    static GrowingVectorMemory<typename Solver::vector_type> vector_memory;
+    static GrowingVectorMemory<Range> vector_memory;
 
-    Vector *v2 = vector_memory.alloc();
+    Range *v2 = vector_memory.alloc();
     op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/ false);
     solver.solve(transpose_operator(op), *v2, u, preconditioner);
     v += *v2;
@@ -696,11 +754,12 @@ inverse_operator(const LinearOperator<typename Solver::vector_type, typename Sol
  *
  * @ingroup LAOperators
  */
-template <typename Range>
-LinearOperator<Range, Range>
+template <typename Range,
+          typename Payload = internal::LinearOperator::EmptyPayload>
+LinearOperator<Range, Range, Payload>
 identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
 {
-  LinearOperator<Range, Range> return_op;
+  LinearOperator<Range, Range, Payload> return_op ((Payload()));
 
   return_op.reinit_range_vector = reinit_vector;
   return_op.reinit_domain_vector = reinit_vector;
@@ -738,9 +797,9 @@ identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
-LinearOperator<Range, Domain>
-null_operator(const LinearOperator<Range, Domain> &op)
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+null_operator(const LinearOperator<Range, Domain, Payload> &op)
 {
   auto return_op = op;
 
@@ -821,6 +880,77 @@ namespace internal
         v.reinit(matrix.n(), omit_zeroing_entries);
       }
     };
+
+
+    /**
+    * A dummy class for LinearOperators that do not require any extensions
+    * to facilitate the operations of the matrix.
+    *
+    * This is the Payload class typically associated with deal.II's native
+    * SparseMatrix. To use Trilinos and PETSc SparseMatrices it is necessary
+    * to initialize a LinearOperator with their associated Payload.
+    *
+    * @author Jean-Paul Pelteret, Matthias Maier, 2016
+    *
+    * @ingroup LAOperators
+    */
+    class EmptyPayload
+    {
+    public:
+
+      /**
+      * Default constructor
+      *
+      * Since this class does not do anything in particular and needs no special
+      * configuration, we have only one generic constructor that can be called
+      * under any conditions.
+      */
+      template <typename... Args>
+      EmptyPayload (const Args &...)
+      { }
+
+
+      /**
+      * Returns a payload configured for transpose operations
+      */
+      EmptyPayload
+      transpose_payload () const
+      {
+        return *this;
+      }
+
+
+      /**
+      * Returns a payload configured for inverse operations
+      */
+      template <typename Solver, typename Preconditioner>
+      EmptyPayload
+      inverse_payload (Solver &, const Preconditioner &) const
+      {
+        return *this;
+      }
+    };
+
+    /**
+    * Operator that returns a payload configured to support the
+    * addition of two LinearOperators
+    */
+    EmptyPayload operator+(const EmptyPayload &,
+                           const EmptyPayload &)
+    {
+      return EmptyPayload();
+    }
+
+    /**
+    * Operator that returns a payload configured to support the
+    * multiplication of two LinearOperators
+    */
+    EmptyPayload operator*(const EmptyPayload &,
+                           const EmptyPayload &)
+    {
+      return EmptyPayload();
+    }
+
   } /* namespace LinearOperator */
 } /* namespace internal */
 
@@ -873,12 +1003,12 @@ namespace
 
   // A helper class to add a reduced matrix interface to a LinearOperator
   // (typically provided by Preconditioner classes)
-  template <typename Range, typename Domain>
+  template <typename Range, typename Domain, typename Payload>
   class MatrixInterfaceWithoutVmultAdd
   {
   public:
     template <typename Matrix>
-    void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+    void operator()(LinearOperator<Range, Domain, Payload> &op, const Matrix &matrix)
     {
       op.vmult = [&matrix](Range &v, const Domain &u)
       {
@@ -939,16 +1069,16 @@ namespace
 
 
   // A helper class to add the full matrix interface to a LinearOperator
-  template <typename Range, typename Domain>
+  template <typename Range, typename Domain, typename Payload>
   class MatrixInterfaceWithVmultAdd
   {
   public:
     template <typename Matrix>
-    void operator()(LinearOperator<Range, Domain> &op, const Matrix &matrix)
+    void operator()(LinearOperator<Range, Domain, Payload> &op, const Matrix &matrix)
     {
       // As above ...
 
-      MatrixInterfaceWithoutVmultAdd<Range, Domain>().operator()(op, matrix);
+      MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>().operator()(op, matrix);
 
       // ... but add native vmult_add and Tvmult_add variants:
 
@@ -1046,11 +1176,13 @@ namespace
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain, typename Matrix>
-LinearOperator<Range, Domain> linear_operator(const Matrix &matrix)
+template <typename Range, typename Domain,
+          typename Payload, typename Matrix>
+LinearOperator<Range, Domain, Payload>
+linear_operator(const Matrix &matrix)
 {
   // implement with the more generic variant below...
-  return linear_operator<Range, Domain, Matrix, Matrix>(matrix, matrix);
+  return linear_operator<Range, Domain, Payload, Matrix, Matrix>(matrix, matrix);
 }
 
 
@@ -1071,12 +1203,14 @@ LinearOperator<Range, Domain> linear_operator(const Matrix &matrix)
  */
 template <typename Range,
           typename Domain,
+          typename Payload,
           typename OperatorExemplar,
           typename Matrix>
-LinearOperator<Range, Domain>
+LinearOperator<Range, Domain, Payload>
 linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
 {
-  LinearOperator<Range, Domain> return_op;
+  // Initialise the payload based on the input exemplar matrix
+  LinearOperator<Range, Domain, Payload> return_op (Payload(operator_exemplar,matrix));
 
   // Always store a reference to matrix and operator_exemplar in the lambda
   // functions. This ensures that a modification of the matrix after the
@@ -1095,15 +1229,17 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix)
 
   typename std::conditional<
   has_vmult_add_and_Tvmult_add<Range, Domain, Matrix>::type::value,
-                               MatrixInterfaceWithVmultAdd<Range, Domain>,
-                               MatrixInterfaceWithoutVmultAdd<Range, Domain>>::type().
+                               MatrixInterfaceWithVmultAdd<Range, Domain, Payload>,
+                               MatrixInterfaceWithoutVmultAdd<Range, Domain, Payload>>::type().
                                operator()(return_op, matrix);
 
   return return_op;
 }
 
+
 //@}
 
+
 DEAL_II_NAMESPACE_CLOSE
 
 #endif // DEAL_II_WITH_CXX11
diff --git a/include/deal.II/lac/packaged_operation.h b/include/deal.II/lac/packaged_operation.h
index 6750437a17..16c3b241f7 100644
--- a/include/deal.II/lac/packaged_operation.h
+++ b/include/deal.II/lac/packaged_operation.h
@@ -28,7 +28,7 @@ DEAL_II_NAMESPACE_OPEN
 
 // Forward declarations:
 template <typename Number> class Vector;
-template <typename Range, typename Domain> class LinearOperator;
+template <typename Range, typename Domain, typename Payload> class LinearOperator;
 template <typename Range = Vector<double> > class PackagedOperation;
 
 
@@ -660,9 +660,9 @@ PackagedOperation<Range> operator*(typename Range::value_type number,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
+template <typename Range, typename Domain, typename Payload>
 PackagedOperation<Range>
-operator*(const LinearOperator<Range, Domain> &op,
+operator*(const LinearOperator<Range, Domain, Payload> &op,
           const Domain &u)
 {
   PackagedOperation<Range> return_comp;
@@ -702,10 +702,10 @@ operator*(const LinearOperator<Range, Domain> &op,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
+template <typename Range, typename Domain, typename Payload>
 PackagedOperation<Domain>
 operator*(const Range &u,
-          const LinearOperator<Range, Domain> &op)
+          const LinearOperator<Range, Domain, Payload> &op)
 {
   PackagedOperation<Range> return_comp;
 
@@ -736,9 +736,9 @@ operator*(const Range &u,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
+template <typename Range, typename Domain, typename Payload>
 PackagedOperation<Range>
-operator*(const LinearOperator<Range, Domain> &op,
+operator*(const LinearOperator<Range, Domain, Payload> &op,
           const PackagedOperation<Domain> &comp)
 {
   PackagedOperation<Range> return_comp;
@@ -786,10 +786,10 @@ operator*(const LinearOperator<Range, Domain> &op,
  *
  * @ingroup LAOperators
  */
-template <typename Range, typename Domain>
+template <typename Range, typename Domain, typename Payload>
 PackagedOperation<Domain>
 operator*(const PackagedOperation<Range> &comp,
-          const LinearOperator<Range, Domain> &op)
+          const LinearOperator<Range, Domain, Payload> &op)
 {
   PackagedOperation<Range> return_comp;
 
diff --git a/include/deal.II/lac/schur_complement.h b/include/deal.II/lac/schur_complement.h
index 1928ee3bb0..17c9d50d07 100644
--- a/include/deal.II/lac/schur_complement.h
+++ b/include/deal.II/lac/schur_complement.h
@@ -224,14 +224,15 @@ DEAL_II_NAMESPACE_OPEN
  * @ingroup LAOperators
  */
 template <typename Range_1, typename Domain_1,
-          typename Range_2, typename Domain_2>
-LinearOperator<Range_2, Domain_2>
-schur_complement(const LinearOperator<Domain_1, Range_1> &A_inv,
-                 const LinearOperator<Range_1, Domain_2> &B,
-                 const LinearOperator<Range_2, Domain_1> &C,
-                 const LinearOperator<Range_2, Domain_2> &D)
+          typename Range_2, typename Domain_2,
+          typename Payload>
+LinearOperator<Range_2, Domain_2, Payload>
+schur_complement(const LinearOperator<Domain_1, Range_1, Payload> &A_inv,
+                 const LinearOperator<Range_1, Domain_2, Payload> &B,
+                 const LinearOperator<Range_2, Domain_1, Payload> &C,
+                 const LinearOperator<Range_2, Domain_2, Payload> &D)
 {
-  LinearOperator<Range_2, Domain_2> return_op;
+  LinearOperator<Range_2, Domain_2, Payload> return_op ((Payload(D)));
 
   return_op.reinit_range_vector = D.reinit_range_vector;
   return_op.reinit_domain_vector = D.reinit_domain_vector;
@@ -361,12 +362,12 @@ schur_complement(const LinearOperator<Domain_1, Range_1> &A_inv,
  * @ingroup LAOperators
  */
 template <typename Range_1, typename Domain_1,
-          typename Range_2>
+          typename Range_2, typename Payload>
 PackagedOperation<Range_2>
-condense_schur_rhs (const LinearOperator<Range_1, Domain_1> &A_inv,
-                    const LinearOperator<Range_2, Domain_1> &C,
-                    const Range_1                           &f,
-                    const Range_2                           &g)
+condense_schur_rhs (const LinearOperator<Range_1, Domain_1, Payload> &A_inv,
+                    const LinearOperator<Range_2, Domain_1, Payload> &C,
+                    const Range_1                                    &f,
+                    const Range_2                                    &g)
 {
   PackagedOperation<Range_2> return_comp;
 
@@ -442,12 +443,12 @@ condense_schur_rhs (const LinearOperator<Range_1, Domain_1> &A_inv,
  * @ingroup LAOperators
  */
 template <typename Range_1, typename Domain_1,
-          typename Domain_2>
+          typename Domain_2, typename Payload>
 PackagedOperation<Domain_1>
-postprocess_schur_solution (const LinearOperator<Range_1, Domain_1> &A_inv,
-                            const LinearOperator<Range_1, Domain_2> &B,
-                            const Domain_2                          &y,
-                            const Range_1                           &f)
+postprocess_schur_solution (const LinearOperator<Range_1, Domain_1, Payload> &A_inv,
+                            const LinearOperator<Range_1, Domain_2, Payload> &B,
+                            const Domain_2                                   &y,
+                            const Range_1                                    &f)
 {
   PackagedOperation<Domain_1> return_comp;
 
diff --git a/tests/lac/linear_operator_01.cc b/tests/lac/linear_operator_01.cc
index 3df7c57d40..f2058ce1fd 100644
--- a/tests/lac/linear_operator_01.cc
+++ b/tests/lac/linear_operator_01.cc
@@ -102,7 +102,8 @@ int main()
 
   // Create to distinct linear operators:
 
-  LinearOperator<LeftVector, RightVector> multiply2;
+  typedef dealii::internal::LinearOperator::EmptyPayload Payload;
+  LinearOperator<LeftVector, RightVector, Payload> multiply2;
   multiply2.vmult = [](LeftVector &v, const RightVector &u)
   {
     v.value = 2 * u.value;
diff --git a/tests/lac/linear_operator_01a.cc b/tests/lac/linear_operator_01a.cc
index b51d8a099e..afb5e1ec1d 100644
--- a/tests/lac/linear_operator_01a.cc
+++ b/tests/lac/linear_operator_01a.cc
@@ -105,7 +105,8 @@ int main()
 
   // Create to distinct linear operators:
 
-  LinearOperator<LeftVector, RightVector> multiply2;
+  typedef dealii::internal::LinearOperator::EmptyPayload Payload;
+  LinearOperator<LeftVector, RightVector, Payload> multiply2;
   multiply2.vmult = [](LeftVector &v, const RightVector &u)
   {
     v.value = 2. * u.value;