lac/linear_operator.h: Add rvalue-reference variant of inverse_operator

diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h
index ef84db6e253c892894817b7c2a9c7e41956aa0d3..f0115ff2578dc3e8ff8713458cd1c76aea8daaf0 100644 (file)
--- a/include/deal.II/lac/linear_operator.h
+++ b/include/deal.II/lac/linear_operator.h
@@ -41,6 +41,8 @@ namespace internal
 template <typename Number>
 class Vector;
 
+class PreconditionIdentity;
+
 template <typename Range  = Vector<double>,
           typename Domain = Range,
           typename Payload =
@@ -71,6 +73,10 @@ template <
 LinearOperator<Range, Domain, Payload>
 null_operator(const LinearOperator<Range, Domain, Payload> &);
 
+template <typename Range, typename Domain, typename Payload>
+LinearOperator<Range, Domain, Payload>
+identity_operator(const LinearOperator<Range, Domain, Payload> &);
+
 
 /**
  * A class to store the abstract concept of a linear operator.
@@ -629,6 +635,7 @@ operator*(const LinearOperator<Range, Intermediate, Payload> & first_op,
     }
 }
 
+
 /**
  * @relatesalso LinearOperator
  *
@@ -725,6 +732,105 @@ inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
   return return_op;
 }
 
+
+/**
+ * @relatesalso LinearOperator
+ *
+ * Variant of above function that takes a LinearOperator @p preconditioner
+ * as preconditioner argument.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Payload,
+          typename Solver,
+          typename Range  = typename Solver::vector_type,
+          typename Domain = Range>
+LinearOperator<Domain, Range, Payload>
+inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
+                 Solver &                                      solver,
+                 const LinearOperator<Range, Domain, Payload> &preconditioner)
+{
+  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](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](Range &       v,
+                                                      const Domain &u) {
+    GrowingVectorMemory<Range> vector_memory;
+
+    typename VectorMemory<Range>::Pointer v2(vector_memory);
+    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
+    solver.solve(op, *v2, u, preconditioner);
+    v += *v2;
+  };
+
+  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](Range &       v,
+                                                       const Domain &u) {
+    GrowingVectorMemory<Range> vector_memory;
+
+    typename VectorMemory<Range>::Pointer v2(vector_memory);
+    op.reinit_range_vector(*v2, /*bool omit_zeroing_entries =*/false);
+    solver.solve(transpose_operator(op), *v2, u, preconditioner);
+    v += *v2;
+  };
+
+  return return_op;
+}
+
+
+/**
+ * @relatesalso LinearOperator
+ *
+ * Variant of above function without a preconditioner argument. In this
+ * case the identity_operator() of the @p op argument is used as a
+ * preconditioner. This is equivalent to using PreconditionIdentity.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Payload,
+          typename Solver,
+          typename Range  = typename Solver::vector_type,
+          typename Domain = Range>
+LinearOperator<Domain, Range, Payload>
+inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
+                 Solver &                                      solver)
+{
+  return inverse_operator(op, solver, identity_operator(op));
+}
+
+
+/**
+ * @relatesalso LinearOperator
+ *
+ * Special overload of above function that takes a PreconditionIdentity
+ * argument.
+ *
+ * @ingroup LAOperators
+ */
+template <typename Payload,
+          typename Solver,
+          typename Range  = typename Solver::vector_type,
+          typename Domain = Range>
+LinearOperator<Domain, Range, Payload>
+inverse_operator(const LinearOperator<Range, Domain, Payload> &op,
+                 Solver &                                      solver,
+                 const PreconditionIdentity &)
+{
+  return inverse_operator(op, solver);
+}
+
 //@}