From: Matthias Maier Date: Tue, 12 Feb 2019 20:34:17 +0000 (-0600) Subject: lac/linear_operator.h: Add rvalue-reference variant of inverse_operator X-Git-Tag: v9.1.0-rc1~332^2~5 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c7a3e6a787f0a1ddc8eb5ae5d43a77df2680eec8;p=dealii.git 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 ef84db6e25..f0115ff257 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -41,6 +41,8 @@ namespace internal template class Vector; +class PreconditionIdentity; + template , typename Domain = Range, typename Payload = @@ -71,6 +73,10 @@ template < LinearOperator null_operator(const LinearOperator &); +template +LinearOperator +identity_operator(const LinearOperator &); + /** * A class to store the abstract concept of a linear operator. @@ -629,6 +635,7 @@ operator*(const LinearOperator & first_op, } } + /** * @relatesalso LinearOperator * @@ -725,6 +732,105 @@ inverse_operator(const LinearOperator &op, return return_op; } + +/** + * @relatesalso LinearOperator + * + * Variant of above function that takes a LinearOperator @p preconditioner + * as preconditioner argument. + * + * @ingroup LAOperators + */ +template +LinearOperator +inverse_operator(const LinearOperator &op, + Solver & solver, + const LinearOperator &preconditioner) +{ + LinearOperator 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 vector_memory; + + typename VectorMemory::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 vector_memory; + + typename VectorMemory::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 +LinearOperator +inverse_operator(const LinearOperator &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 +LinearOperator +inverse_operator(const LinearOperator &op, + Solver & solver, + const PreconditionIdentity &) +{ + return inverse_operator(op, solver); +} + //@}