From: Jean-Paul Pelteret 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. +
+(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 + class EmptyBlockPayload; + } +} + template class BlockVector; template , - typename Domain = Range> + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > class BlockLinearOperator; template , typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>, typename BlockMatrixType> -BlockLinearOperator +BlockLinearOperator block_operator(const BlockMatrixType &matrix); template , - typename Domain = Range> -BlockLinearOperator -block_operator(const std::array, n>, m> &); + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > +BlockLinearOperator +block_operator(const std::array, n>, m> &); template , - typename Domain = Range> -BlockLinearOperator -block_diagonal_operator(const std::array, m> &); + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > +BlockLinearOperator +block_diagonal_operator(const std::array, m> &); template , - typename Domain = Range> -BlockLinearOperator -block_diagonal_operator(const LinearOperator &op); + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > +BlockLinearOperator +block_diagonal_operator(const LinearOperator &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 Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<>, typename BlockMatrixType> -BlockLinearOperator +BlockLinearOperator block_diagonal_operator(const BlockMatrixType &block_matrix); template , - typename Domain = Range> -LinearOperator -block_forward_substitution(const BlockLinearOperator &, - const BlockLinearOperator &); + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > +LinearOperator +block_forward_substitution(const BlockLinearOperator &, + const BlockLinearOperator &); template , - typename Domain = Range> -LinearOperator -block_back_substitution(const BlockLinearOperator &, - const BlockLinearOperator &); + typename Domain = Range, + typename BlockPayload = internal::BlockLinearOperator::EmptyBlockPayload<> > +LinearOperator +block_back_substitution(const BlockLinearOperator &, + const BlockLinearOperator &); // end of workaround @@ -125,12 +142,12 @@ block_back_substitution(const BlockLinearOperator &, * * @ingroup LAOperators */ -template -class BlockLinearOperator : public LinearOperator +template +class BlockLinearOperator : public LinearOperator { public: - typedef LinearOperator BlockType; + typedef LinearOperator 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() + BlockLinearOperator(const BlockPayload &payload) + : LinearOperator(typename BlockPayload::BlockType(payload,payload)) { - n_block_rows = []() -> unsigned int { Assert(false, ExcMessage("Uninitialized BlockLinearOperator::n_block_rows called")); @@ -165,7 +181,7 @@ public: /** * Default copy constructor. */ - BlockLinearOperator(const BlockLinearOperator &) = + BlockLinearOperator(const BlockLinearOperator &) = default; /** @@ -176,7 +192,7 @@ public: template BlockLinearOperator(const Op &op) { - *this = block_operator(op); + *this = block_operator(op); } /** @@ -187,7 +203,7 @@ public: template BlockLinearOperator(const std::array, m> &ops) { - *this = block_operator(ops); + *this = block_operator(ops); } /** @@ -198,23 +214,23 @@ public: template BlockLinearOperator(const std::array &ops) { - *this = block_diagonal_operator(ops); + *this = block_diagonal_operator(ops); } /** * Default copy assignment operator. */ - BlockLinearOperator & - operator=(const BlockLinearOperator &) = default; + BlockLinearOperator & + operator=(const BlockLinearOperator &) = default; /** * Templated copy assignment operator for an object @p op for which the * conversion function block_operator is defined. */ template - BlockLinearOperator &operator=(const Op &op) + BlockLinearOperator &operator=(const Op &op) { - *this = block_operator(op); + *this = block_operator(op); return *this; } @@ -224,10 +240,10 @@ public: * specialization. */ template - BlockLinearOperator & + BlockLinearOperator & operator=(const std::array, m> &ops) { - *this = block_operator(ops); + *this = block_operator(ops); return *this; } @@ -237,10 +253,10 @@ public: * operator calls the corresponding block_operator() specialization. */ template - BlockLinearOperator & + BlockLinearOperator & operator=(const std::array &ops) { - *this = block_diagonal_operator(ops); + *this = block_diagonal_operator(ops); return *this; } @@ -272,10 +288,10 @@ namespace internal { // Populate the LinearOperator interfaces with the help of the // BlockLinearOperator functions - template + template inline void populate_linear_operator_functions( - dealii::BlockLinearOperator &op) + dealii::BlockLinearOperator &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 + 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 + EmptyBlockPayload (const Args &...) + { } + }; + } /*namespace BlockLinearOperator*/ } /*namespace internal*/ @@ -382,13 +434,14 @@ namespace internal */ template -BlockLinearOperator +BlockLinearOperator block_operator(const BlockMatrixType &block_matrix) { - typedef typename BlockLinearOperator::BlockType BlockType; + typedef typename BlockLinearOperator::BlockType BlockType; - BlockLinearOperator return_op; + BlockLinearOperator 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 -BlockLinearOperator -block_operator(const std::array, n>, m> &ops) +template +BlockLinearOperator +block_operator(const std::array, n>, m> &ops) { static_assert(m > 0 && n > 0, "a blocked LinearOperator must consist of at least one block"); - typedef typename BlockLinearOperator::BlockType BlockType; + typedef typename BlockLinearOperator::BlockType BlockType; - BlockLinearOperator return_op; + BlockLinearOperator 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 -BlockLinearOperator +BlockLinearOperator block_diagonal_operator(const BlockMatrixType &block_matrix) { - typedef typename BlockLinearOperator::BlockType BlockType; + typedef typename BlockLinearOperator::BlockType BlockType; - BlockLinearOperator return_op; + BlockLinearOperator 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 -BlockLinearOperator -block_diagonal_operator(const std::array, m> &ops) +template +BlockLinearOperator +block_diagonal_operator(const std::array, m> &ops) { static_assert(m > 0, "a blockdiagonal LinearOperator must consist of at least one block"); - typedef typename BlockLinearOperator::BlockType BlockType; + typedef typename BlockLinearOperator::BlockType BlockType; std::array, m> new_ops; @@ -597,15 +651,15 @@ block_diagonal_operator(const std::array -BlockLinearOperator -block_diagonal_operator(const LinearOperator &op) +template +BlockLinearOperator +block_diagonal_operator(const LinearOperator &op) { static_assert(m > 0, "a blockdiagonal LinearOperator must consist of at least " "one block"); - typedef typename BlockLinearOperator::BlockType BlockType; + typedef typename BlockLinearOperator::BlockType BlockType; std::array new_ops; new_ops.fill(op); @@ -655,12 +709,12 @@ block_diagonal_operator(const LinearOperator -LinearOperator -block_forward_substitution(const BlockLinearOperator &block_operator, - const BlockLinearOperator &diagonal_inverse) +template +LinearOperator +block_forward_substitution(const BlockLinearOperator &block_operator, + const BlockLinearOperator &diagonal_inverse) { - LinearOperator return_op; + LinearOperator 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 &block_opera * * @ingroup LAOperators */ -template -LinearOperator -block_back_substitution(const BlockLinearOperator &block_operator, - const BlockLinearOperator &diagonal_inverse) +template +LinearOperator +block_back_substitution(const BlockLinearOperator &block_operator, + const BlockLinearOperator &diagonal_inverse) { - LinearOperator return_op; + LinearOperator 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 class Vector; -template , typename Domain = Range> +template , + typename Domain = Range, + typename Payload = internal::LinearOperator::EmptyPayload> class LinearOperator; template , typename Domain = Range, + typename Payload = internal::LinearOperator::EmptyPayload, typename OperatorExemplar, typename Matrix> -LinearOperator linear_operator (const OperatorExemplar &, - const Matrix &); +LinearOperator linear_operator (const OperatorExemplar &, + const Matrix &); template , typename Domain = Range, + typename Payload = internal::LinearOperator::EmptyPayload, typename Matrix> -LinearOperator linear_operator (const Matrix &); +LinearOperator linear_operator (const Matrix &); template , - typename Domain = Range> -LinearOperator -null_operator(const LinearOperator &); + typename Domain = Range, + typename Payload = internal::LinearOperator::EmptyPayload> +LinearOperator +null_operator(const LinearOperator &); /** @@ -104,25 +117,58 @@ null_operator(const LinearOperator &); * 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 DEAL_II_WITH_CXX11 is enabled during cmake * configure. * - * @author Luca Heltai, Matthias Maier, 2015 + * @author Luca Heltai, Matthias Maier, 2015; Jean-Paul Pelteret, 2016 * * @ingroup LAOperators */ -template class LinearOperator +template +class LinearOperator + : public Payload { public: /** - * Create an empty LinearOperator object. All std::function - * 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 &) = default; + LinearOperator (const LinearOperator &) = default; /** * Templated copy constructor that creates a LinearOperator object from an @@ -174,27 +220,27 @@ public: * linear_operator is defined. */ template, Op>::value>::type> + typename = typename std::enable_if, Op>::value>::type> LinearOperator (const Op &op) { - *this = linear_operator(op); + *this = linear_operator(op); } /** * Default copy assignment operator. */ - LinearOperator & - operator=(const LinearOperator &) = default; + LinearOperator & + operator=(const LinearOperator &) = default; /** * Templated copy assignment operator for an object @p op for which the * conversion function linear_operator is defined. */ template , Op>::value>::type> - LinearOperator &operator=(const Op &op) + typename = typename std::enable_if, Op>::value>::type> + LinearOperator &operator=(const Op &op) { - *this = linear_operator(op); + *this = linear_operator(op); return *this; } @@ -249,8 +295,8 @@ public: * Addition with a LinearOperator @p second_op with the same @p Domain and * @p Range. */ - LinearOperator & - operator+=(const LinearOperator &second_op) + LinearOperator & + operator+=(const LinearOperator &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 & - operator-=(const LinearOperator &second_op) + LinearOperator & + operator-=(const LinearOperator &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 & - operator*=(const LinearOperator &second_op) + LinearOperator & + operator*=(const LinearOperator &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 + LinearOperator operator*=(typename Domain::value_type number) { *this = *this * number; @@ -314,10 +360,10 @@ public: * * @ingroup LAOperators */ -template -LinearOperator -operator+(const LinearOperator &first_op, - const LinearOperator &second_op) +template +LinearOperator +operator+(const LinearOperator &first_op, + const LinearOperator &second_op) { if (first_op.is_null_operator) { @@ -329,7 +375,9 @@ operator+(const LinearOperator &first_op, } else { - LinearOperator return_op; + LinearOperator return_op ( + static_cast(first_op) + static_cast(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 &first_op, * * @ingroup LAOperators */ -template -LinearOperator -operator-(const LinearOperator &first_op, - const LinearOperator &second_op) +template +LinearOperator +operator-(const LinearOperator &first_op, + const LinearOperator &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 &first_op, * * @ingroup LAOperators */ -template -LinearOperator +template +LinearOperator operator*(typename Range::value_type number, - const LinearOperator &op) + const LinearOperator &op) { static_assert( std::is_convertible::value, @@ -432,7 +480,7 @@ operator*(typename Range::value_type number, } else { - LinearOperator return_op = op; + LinearOperator 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 -LinearOperator -operator*(const LinearOperator &op, +template +LinearOperator +operator*(const LinearOperator &op, typename Domain::value_type number) { static_assert( @@ -512,21 +560,23 @@ operator*(const LinearOperator &op, * * @ingroup LAOperators */ -template -LinearOperator -operator*(const LinearOperator &first_op, - const LinearOperator &second_op) +template +LinearOperator +operator*(const LinearOperator &first_op, + const LinearOperator &second_op) { if (first_op.is_null_operator || second_op.is_null_operator) { - LinearOperator return_op; + LinearOperator 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 return_op; + LinearOperator return_op ( + static_cast(first_op) * static_cast(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 &first_op, * * Return the transpose linear operations of @p op. * + * @author Matthias Maier, 2015 + * * @ingroup LAOperators */ -template -LinearOperator -transpose_operator(const LinearOperator &op) +template +LinearOperator +transpose_operator(const LinearOperator &op) { - LinearOperator return_op; + LinearOperator 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 &op) return return_op; } + /** * @relates LinearOperator * @@ -622,52 +677,55 @@ transpose_operator(const LinearOperator &op) * of the @p solver object will be modified upon invocation of * vmult or Tvmult. * + * @author Luca Heltai, Matthias Maier, Jean-Paul Pelteret, 2015 + * * @ingroup LAOperators */ -template -LinearOperator -inverse_operator(const LinearOperator &op, +template +LinearOperator +inverse_operator(const LinearOperator &op, Solver &solver, const Preconditioner &preconditioner) { - typedef typename Solver::vector_type Vector; - - LinearOperator return_op; + 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](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 vector_memory; + static GrowingVectorMemory 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 vector_memory; + static GrowingVectorMemory 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 -LinearOperator +template +LinearOperator identity_operator(const std::function &reinit_vector) { - LinearOperator return_op; + LinearOperator 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 &reinit_vector) * * @ingroup LAOperators */ -template -LinearOperator -null_operator(const LinearOperator &op) +template +LinearOperator +null_operator(const LinearOperator &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 + EmptyPayload (const Args &...) + { } + + + /** + * Returns a payload configured for transpose operations + */ + EmptyPayload + transpose_payload () const + { + return *this; + } + + + /** + * Returns a payload configured for inverse operations + */ + template + 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 + template class MatrixInterfaceWithoutVmultAdd { public: template - void operator()(LinearOperator &op, const Matrix &matrix) + void operator()(LinearOperator &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 + template class MatrixInterfaceWithVmultAdd { public: template - void operator()(LinearOperator &op, const Matrix &matrix) + void operator()(LinearOperator &op, const Matrix &matrix) { // As above ... - MatrixInterfaceWithoutVmultAdd().operator()(op, matrix); + MatrixInterfaceWithoutVmultAdd().operator()(op, matrix); // ... but add native vmult_add and Tvmult_add variants: @@ -1046,11 +1176,13 @@ namespace * * @ingroup LAOperators */ -template -LinearOperator linear_operator(const Matrix &matrix) +template +LinearOperator +linear_operator(const Matrix &matrix) { // implement with the more generic variant below... - return linear_operator(matrix, matrix); + return linear_operator(matrix, matrix); } @@ -1071,12 +1203,14 @@ LinearOperator linear_operator(const Matrix &matrix) */ template -LinearOperator +LinearOperator linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix) { - LinearOperator return_op; + // Initialise the payload based on the input exemplar matrix + LinearOperator 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::type::value, - MatrixInterfaceWithVmultAdd, - MatrixInterfaceWithoutVmultAdd>::type(). + MatrixInterfaceWithVmultAdd, + MatrixInterfaceWithoutVmultAdd>::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 class Vector; -template class LinearOperator; +template class LinearOperator; template > class PackagedOperation; @@ -660,9 +660,9 @@ PackagedOperation operator*(typename Range::value_type number, * * @ingroup LAOperators */ -template +template PackagedOperation -operator*(const LinearOperator &op, +operator*(const LinearOperator &op, const Domain &u) { PackagedOperation return_comp; @@ -702,10 +702,10 @@ operator*(const LinearOperator &op, * * @ingroup LAOperators */ -template +template PackagedOperation operator*(const Range &u, - const LinearOperator &op) + const LinearOperator &op) { PackagedOperation return_comp; @@ -736,9 +736,9 @@ operator*(const Range &u, * * @ingroup LAOperators */ -template +template PackagedOperation -operator*(const LinearOperator &op, +operator*(const LinearOperator &op, const PackagedOperation &comp) { PackagedOperation return_comp; @@ -786,10 +786,10 @@ operator*(const LinearOperator &op, * * @ingroup LAOperators */ -template +template PackagedOperation operator*(const PackagedOperation &comp, - const LinearOperator &op) + const LinearOperator &op) { PackagedOperation 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 -LinearOperator -schur_complement(const LinearOperator &A_inv, - const LinearOperator &B, - const LinearOperator &C, - const LinearOperator &D) + typename Range_2, typename Domain_2, + typename Payload> +LinearOperator +schur_complement(const LinearOperator &A_inv, + const LinearOperator &B, + const LinearOperator &C, + const LinearOperator &D) { - LinearOperator return_op; + LinearOperator 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 &A_inv, * @ingroup LAOperators */ template + typename Range_2, typename Payload> PackagedOperation -condense_schur_rhs (const LinearOperator &A_inv, - const LinearOperator &C, - const Range_1 &f, - const Range_2 &g) +condense_schur_rhs (const LinearOperator &A_inv, + const LinearOperator &C, + const Range_1 &f, + const Range_2 &g) { PackagedOperation return_comp; @@ -442,12 +443,12 @@ condense_schur_rhs (const LinearOperator &A_inv, * @ingroup LAOperators */ template + typename Domain_2, typename Payload> PackagedOperation -postprocess_schur_solution (const LinearOperator &A_inv, - const LinearOperator &B, - const Domain_2 &y, - const Range_1 &f) +postprocess_schur_solution (const LinearOperator &A_inv, + const LinearOperator &B, + const Domain_2 &y, + const Range_1 &f) { PackagedOperation 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 multiply2; + typedef dealii::internal::LinearOperator::EmptyPayload Payload; + LinearOperator 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 multiply2; + typedef dealii::internal::LinearOperator::EmptyPayload Payload; + LinearOperator multiply2; multiply2.vmult = [](LeftVector &v, const RightVector &u) { v.value = 2. * u.value;