From 8cea49640f44929bf8b744d2e88284123006108f Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Mon, 27 Apr 2015 12:02:29 +0200 Subject: [PATCH] Restructure documentation --- include/deal.II/lac/computation.h | 24 +- include/deal.II/lac/linear_operator.h | 804 +++++++++++++------------- 2 files changed, 428 insertions(+), 400 deletions(-) diff --git a/include/deal.II/lac/computation.h b/include/deal.II/lac/computation.h index bcadf0fd75..a583ef267d 100644 --- a/include/deal.II/lac/computation.h +++ b/include/deal.II/lac/computation.h @@ -265,7 +265,7 @@ public: //@{ /** - * \relates Computation + * @relates Computation * * Addition of two Computation objects @p first_comp and @p second_comp given by * vector space addition of the corresponding results. @@ -300,7 +300,7 @@ operator+(const Computation &first_comp, } /** - * \relates Computation + * @relates Computation * * Subtraction of two Computation objects @p first_comp and @p second_comp * given by vector space addition of the corresponding results. @@ -338,7 +338,7 @@ operator-(const Computation &first_comp, } /** - * \relates Computation + * @relates Computation * * Scalar multiplication of a Computation objects @p comp with a scalar @p number * given by a scaling Computation result with @p number. @@ -371,7 +371,7 @@ operator*(const Computation &comp, } /** - * \relates Computation + * @relates Computation * * Scalar multiplication of a Computation objects @p comp with a scalar @p number * given by a scaling Computation result with @p number. @@ -387,7 +387,7 @@ operator*(typename Range::value_type number, } /** - * \relates Computation + * @relates Computation * * Add a constant @p offset (of the @p Range space) to the result of a * Computation. @@ -402,7 +402,7 @@ Computation operator+(const Computation &comp, } /** - * \relates Computation + * @relates Computation * * Add a constant @p offset (of the @p Range space) to the result of a * Computation. @@ -417,7 +417,7 @@ Computation operator+(const Range &offset, } /** - * \relates Computation + * @relates Computation * * Subtract a constant @p offset (of the @p Range space) from the result of a * Computation. @@ -433,7 +433,7 @@ Computation operator-(const Computation &comp, /** - * \relates Computation + * @relates Computation * * Subtract a computational result from a constant @p offset (of the @p * Range space). The result is a Computation object that applies this @@ -570,7 +570,7 @@ Computation operator-(const Range &u, const Range &v) /** - * \relates Computation + * @relates Computation * * Create a Computation object from a LinearOperator and a reference to a * vector @p u of the Domain space. The object stores the computation @@ -612,7 +612,7 @@ operator*(const LinearOperator &op, /** - * \relates Computation + * @relates Computation * * Create a Computation object from a LinearOperator and a reference to a * vector @p u of the Range space. The object stores the computation @@ -654,7 +654,7 @@ operator*(const Range &u, /** - * \relates Computation + * @relates Computation * * Composition of a Computation object with a LinearOperator. * The object stores the computation $\text{op} \,comp$ (in matrix notation). @@ -704,7 +704,7 @@ operator*(const LinearOperator &op, /** - * \relates Computation + * @relates Computation * * Composition of a Computation object with a LinearOperator. * The object stores the computation $\text{op}^T \,comp$ (in matrix notation). diff --git a/include/deal.II/lac/linear_operator.h b/include/deal.II/lac/linear_operator.h index c238222d42..71915eb777 100644 --- a/include/deal.II/lac/linear_operator.h +++ b/include/deal.II/lac/linear_operator.h @@ -225,6 +225,10 @@ public: */ std::function reinit_domain_vector; + /** + * @name In-place vector space operations + */ + //@{ /** * Addition with a LinearOperator @p second_op with the same @p Domain @@ -270,12 +274,17 @@ public: return *this; } + //@} }; +/** + * @name Vector space operations + */ +//@{ /** - * \relates LinearOperator + * @relates LinearOperator * * 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)$ @@ -324,7 +333,7 @@ operator+(const LinearOperator &first_op, /** - * \relates LinearOperator + * @relates LinearOperator * * 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)$ @@ -342,76 +351,7 @@ operator-(const LinearOperator &first_op, /** - * \relates LinearOperator - * - * Composition 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))$ - * - * @ingroup LAOperators - */ -template -LinearOperator -operator*(const LinearOperator &first_op, - const LinearOperator &second_op) -{ - LinearOperator 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 - - return_op.vmult = [first_op, second_op](Range &v, const Domain &u) - { - static GrowingVectorMemory vector_memory; - - Intermediate *i = vector_memory.alloc(); - second_op.reinit_range_vector(*i, /*bool fast =*/ true); - second_op.vmult(*i, u); - first_op.vmult(v, *i); - vector_memory.free(i); - }; - - return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) - { - static GrowingVectorMemory vector_memory; - - Intermediate *i = vector_memory.alloc(); - second_op.reinit_range_vector(*i, /*bool fast =*/ true); - second_op.vmult(*i, u); - first_op.vmult_add(v, *i); - vector_memory.free(i); - }; - - return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) - { - static GrowingVectorMemory vector_memory; - - Intermediate *i = vector_memory.alloc(); - first_op.reinit_domain_vector(*i, /*bool fast =*/ true); - first_op.Tvmult(*i, u); - second_op.Tvmult(v, *i); - vector_memory.free(i); - }; - - return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) - { - static GrowingVectorMemory vector_memory; - - Intermediate *i = vector_memory.alloc(); - first_op.reinit_domain_vector(*i, /*bool fast =*/ true); - first_op.Tvmult(*i, u); - second_op.Tvmult_add(v, *i); - vector_memory.free(i); - }; - - return return_op; -} - - -/** - * \relates LinearOperator + * @relates LinearOperator * * Scalar multiplication of a ScalarOperator object @p op with @p * number from the left. @@ -472,7 +412,7 @@ operator*(typename Range::value_type number, /** - * \relates LinearOperator + * @relates LinearOperator * * Scalar multiplication of a ScalarOperator object from the right. * @@ -498,80 +438,108 @@ operator*(const LinearOperator &op, return number * op; } +//@} + /** - * \relates LinearOperator - * - * Returns the transpose linear operations of @ref op. - * - * @ingroup LAOperators + * @name Composition and manipulation of a LinearOperator */ -template -LinearOperator -transpose_operator(const LinearOperator &op) -{ - LinearOperator 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; -} - +//@{ /** - * \relates LinearOperator - * - * Returns a LinearOperator that is the identity of the vector space - * @p Range. + * @relates LinearOperator * - * The function takes an std::function object @ref - * reinit_vector as an argument to initialize the - * reinit_range_vector and reinit_domain_vector - * objects of the LinearOperator object. + * Composition 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))$ * * @ingroup LAOperators */ -template -LinearOperator -identity_operator(const std::function &reinit_vector) +template +LinearOperator +operator*(const LinearOperator &first_op, + const LinearOperator &second_op) { - LinearOperator return_op; + LinearOperator return_op; - return_op.reinit_range_vector = reinit_vector; - return_op.reinit_domain_vector = reinit_vector; + return_op.reinit_domain_vector = second_op.reinit_domain_vector; + return_op.reinit_range_vector = first_op.reinit_range_vector; - return_op.vmult = [](Range &v, const Range &u) + // 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) { - v = u; + static GrowingVectorMemory vector_memory; + + Intermediate *i = vector_memory.alloc(); + second_op.reinit_range_vector(*i, /*bool fast =*/ true); + second_op.vmult(*i, u); + first_op.vmult(v, *i); + vector_memory.free(i); }; - return_op.vmult_add = [](Range &v, const Range &u) + return_op.vmult_add = [first_op, second_op](Range &v, const Domain &u) { - v += u; + static GrowingVectorMemory vector_memory; + + Intermediate *i = vector_memory.alloc(); + second_op.reinit_range_vector(*i, /*bool fast =*/ true); + second_op.vmult(*i, u); + first_op.vmult_add(v, *i); + vector_memory.free(i); }; - return_op.Tvmult = [](Range &v, const Range &u) + return_op.Tvmult = [first_op, second_op](Domain &v, const Range &u) { - v = u; + static GrowingVectorMemory vector_memory; + + Intermediate *i = vector_memory.alloc(); + first_op.reinit_domain_vector(*i, /*bool fast =*/ true); + first_op.Tvmult(*i, u); + second_op.Tvmult(v, *i); + vector_memory.free(i); }; - return_op.Tvmult_add = [](Range &v, const Range &u) + return_op.Tvmult_add = [first_op, second_op](Domain &v, const Range &u) { - v += u; + static GrowingVectorMemory vector_memory; + + Intermediate *i = vector_memory.alloc(); + first_op.reinit_domain_vector(*i, /*bool fast =*/ true); + first_op.Tvmult(*i, u); + second_op.Tvmult_add(v, *i); + vector_memory.free(i); }; return return_op; } +/** + * @relates LinearOperator + * + * Returns the transpose linear operations of @ref op. + * + * @ingroup LAOperators + */ +template +LinearOperator +transpose_operator(const LinearOperator &op) +{ + LinearOperator 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; +} /** - * \relates LinearOperator + * @relates LinearOperator * * Returns an object representing the inverse of the LinearOperator @p op. * @@ -639,297 +607,60 @@ inverse_operator(const LinearOperator) - * @code - * op_a00 | op_a01 - * | - * --------------- - * | - * op_a10 | op_a11 - * @endcode - * The coresponding block_operator invocation takes the form - * @code - * block_operator<2, 2, BlockVector>({op_a00, op_a01, op_a10, op_a11}); - * @endcode + * The function takes an std::function object @ref + * reinit_vector as an argument to initialize the + * reinit_range_vector and reinit_domain_vector + * objects of the LinearOperator object. * * @ingroup LAOperators */ -template , - typename Domain = Range> -LinearOperator -block_operator(const std::array, n>, m> &ops) +template +LinearOperator +identity_operator(const std::function &reinit_vector) { - static_assert(m > 0 && n > 0, - "a blocked LinearOperator must consist of at least one block"); + LinearOperator return_op; - LinearOperator return_op; + return_op.reinit_range_vector = reinit_vector; + return_op.reinit_domain_vector = reinit_vector; - return_op.reinit_range_vector = [ops](Range &v, bool fast) + return_op.vmult = [](Range &v, const Range &u) { - // Reinitialize the block vector to m blocks: - v.reinit(m); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < m; ++i) - ops[i][0].reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); + v = u; }; - return_op.reinit_domain_vector = [ops](Domain &v, bool fast) + return_op.vmult_add = [](Range &v, const Range &u) { - // Reinitialize the block vector to n blocks: - v.reinit(n); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < n; ++i) - ops[0][i].reinit_domain_vector(v.block(i), fast); + v += u; + }; - v.collect_sizes(); + return_op.Tvmult = [](Range &v, const Range &u) + { + v = u; }; - return_op.vmult = [ops](Range &v, const Domain &u) + return_op.Tvmult_add = [](Range &v, const Range &u) { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); - - for (unsigned int i = 0; i < m; ++i) - { - ops[i][0].vmult(v.block(i), u.block(0)); - for (unsigned int j = 1; j < n; ++j) - ops[i][j].vmult_add(v.block(i), u.block(j)); - } - }; - - return_op.vmult_add = [ops](Range &v, const Domain &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); - - for (unsigned int i = 0; i < m; ++i) - for (unsigned int j = 0; j < n; ++j) - ops[i][j].vmult_add(v.block(i), u.block(j)); - }; - - return_op.Tvmult = [ops](Domain &v, const Range &u) - { - Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < n; ++i) - { - ops[0][i].Tvmult(v.block(i), u.block(0)); - for (unsigned int j = 1; j < m; ++j) - ops[j][i].Tvmult_add(v.block(i), u.block(j)); - } - }; - - return_op.Tvmult_add = [ops](Domain &v, const Range &u) - { - Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < n; ++i) - for (unsigned int j = 0; j < m; ++j) - ops[j][i].Tvmult_add(v.block(i), u.block(j)); - }; - - return return_op; -} - - -/** - * \relates LinearOperator - * - * A variant of above function that builds up a block diagonal linear - * operator from an array @p ops of diagonal elements (off-diagonal blocks - * are assumed to be 0). - * - * The list @p ops is best passed as an initializer list. Consider for - * example a linear operator block (acting on Vector) - * diag(op_a0, op_a1, ..., op_am). The coresponding - * block_operator invocation takes the form - * @code - * block_diagonal_operator>({op_00, op_a1, ..., op_am}); - * @endcode - * - * @ingroup LAOperators - */ -template , - typename Domain = Range> -LinearOperator -block_diagonal_operator(const std::array, m> &ops) -{ - static_assert(m > 0, - "a blockdiagonal LinearOperator must consist of at least one block"); - - LinearOperator return_op; - - return_op.reinit_range_vector = [ops](Range &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(m); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < m; ++i) - ops[i].reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.reinit_domain_vector = [ops](Domain &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(m); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < m; ++i) - ops[i].reinit_domain_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.vmult = [ops](Range &v, const Domain &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - ops[i].vmult(v.block(i), u.block(i)); - }; - - return_op.vmult_add = [ops](Range &v, const Domain &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - ops[i].vmult_add(v.block(i), u.block(i)); - }; - - return_op.Tvmult = [ops](Domain &v, const Range &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - ops[i].Tvmult(v.block(i), u.block(i)); - }; - - return_op.Tvmult_add = [ops](Domain &v, const Range &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - ops[i].Tvmult_add(v.block(i), u.block(i)); - }; + v += u; + }; return return_op; } -/** - * \relates LinearOperator - * - * A variant of above function that only takes a single LinearOperator - * argument @p op and creates a blockdiagonal linear operator with @p m - * copies of it. - * - * @ingroup LAOperators - */ -template , - typename Domain = Range> -LinearOperator -block_diagonal_operator(const LinearOperator &op) -{ - - static_assert(m > 0, - "a blockdiagonal LinearOperator must consist of at least " - "one block"); - - LinearOperator return_op; - - return_op.reinit_range_vector = [op](Range &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(m); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < m; ++i) - op.reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.reinit_domain_vector = [op](Domain &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(m); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < m; ++i) - op.reinit_domain_vector(v.block(i), fast); - - v.collect_sizes(); - }; - - return_op.vmult = [op](Range &v, const Domain &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - op.vmult(v.block(i), u.block(i)); - }; - - return_op.vmult_add = [op](Range &v, const Domain &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - op.vmult_add(v.block(i), u.block(i)); - }; - - return_op.Tvmult = [op](Domain &v, const Range &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - op.Tvmult(v.block(i), u.block(i)); - }; - - return_op.Tvmult_add = [op](Domain &v, const Range &u) - { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - - for (unsigned int i = 0; i < m; ++i) - op.Tvmult_add(v.block(i), u.block(i)); - }; - - return return_op; -} - - - namespace internal { namespace LinearOperator @@ -1100,7 +831,7 @@ namespace /** - * \relates LinearOperator + * @relates LinearOperator * * A function that encapsulates generic @p matrix objects that act on a * compatible Vector type into a LinearOperator. The LinearOperator object @@ -1166,7 +897,7 @@ LinearOperator linear_operator(const Matrix &matrix) /** - * \relates LinearOperator + * @relates LinearOperator * * Variant of above function that takes an operator object @p * operator_exemplar as an additional reference. This object is used to @@ -1215,6 +946,303 @@ linear_operator(const OperatorExemplar &operator_exemplar, const Matrix &matrix) return return_op; } +//@} + + +/** + * @name Creation of LinearOperator block structures + */ +//@{ + +/** + * @relates LinearOperator + * + * A function that encapsulates a given collection @p ops of + * LinearOperators into a block structure. Hereby, it is assumed that Range + * and Domain are blockvectors, i.e., derived from @ref BlockVectorBase. + * The individual linear operators in @p ops must act on a the underlying + * vector type of the block vectors, i.e., on Domain::BlockType yielding a + * result in Range::BlockType. + * + * The list @p ops is best passed as an initializer list. Consider for + * example a linear operator block (acting on Vector) + * @code + * op_a00 | op_a01 + * | + * --------------- + * | + * op_a10 | op_a11 + * @endcode + * The coresponding block_operator invocation takes the form + * @code + * block_operator<2, 2, BlockVector>({op_a00, op_a01, op_a10, op_a11}); + * @endcode + * + * @ingroup LAOperators + */ +template , + typename Domain = Range> +LinearOperator +block_operator(const std::array, n>, m> &ops) +{ + static_assert(m > 0 && n > 0, + "a blocked LinearOperator must consist of at least one block"); + + LinearOperator return_op; + + return_op.reinit_range_vector = [ops](Range &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < m; ++i) + ops[i][0].reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [ops](Domain &v, bool fast) + { + // Reinitialize the block vector to n blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < n; ++i) + ops[0][i].reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [ops](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); + + for (unsigned int i = 0; i < m; ++i) + { + ops[i][0].vmult(v.block(i), u.block(0)); + for (unsigned int j = 1; j < n; ++j) + ops[i][j].vmult_add(v.block(i), u.block(j)); + } + }; + + return_op.vmult_add = [ops](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); + + for (unsigned int i = 0; i < m; ++i) + for (unsigned int j = 0; j < n; ++j) + ops[i][j].vmult_add(v.block(i), u.block(j)); + }; + + return_op.Tvmult = [ops](Domain &v, const Range &u) + { + Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < n; ++i) + { + ops[0][i].Tvmult(v.block(i), u.block(0)); + for (unsigned int j = 1; j < m; ++j) + ops[j][i].Tvmult_add(v.block(i), u.block(j)); + } + }; + + return_op.Tvmult_add = [ops](Domain &v, const Range &u) + { + Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < n; ++i) + for (unsigned int j = 0; j < m; ++j) + ops[j][i].Tvmult_add(v.block(i), u.block(j)); + }; + + return return_op; +} + + +/** + * @relates LinearOperator + * + * A variant of above function that builds up a block diagonal linear + * operator from an array @p ops of diagonal elements (off-diagonal blocks + * are assumed to be 0). + * + * The list @p ops is best passed as an initializer list. Consider for + * example a linear operator block (acting on Vector) + * diag(op_a0, op_a1, ..., op_am). The coresponding + * block_operator invocation takes the form + * @code + * block_diagonal_operator>({op_00, op_a1, ..., op_am}); + * @endcode + * + * @ingroup LAOperators + */ +template , + typename Domain = Range> +LinearOperator +block_diagonal_operator(const std::array, m> &ops) +{ + static_assert(m > 0, + "a blockdiagonal LinearOperator must consist of at least one block"); + + LinearOperator return_op; + + return_op.reinit_range_vector = [ops](Range &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < m; ++i) + ops[i].reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [ops](Domain &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < m; ++i) + ops[i].reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [ops](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + ops[i].vmult(v.block(i), u.block(i)); + }; + + return_op.vmult_add = [ops](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + ops[i].vmult_add(v.block(i), u.block(i)); + }; + + return_op.Tvmult = [ops](Domain &v, const Range &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + ops[i].Tvmult(v.block(i), u.block(i)); + }; + + return_op.Tvmult_add = [ops](Domain &v, const Range &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + ops[i].Tvmult_add(v.block(i), u.block(i)); + }; + + return return_op; +} + + +/** + * @relates LinearOperator + * + * A variant of above function that only takes a single LinearOperator + * argument @p op and creates a blockdiagonal linear operator with @p m + * copies of it. + * + * @ingroup LAOperators + */ +template , + typename Domain = Range> +LinearOperator +block_diagonal_operator(const LinearOperator &op) +{ + + static_assert(m > 0, + "a blockdiagonal LinearOperator must consist of at least " + "one block"); + + LinearOperator return_op; + + return_op.reinit_range_vector = [op](Range &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_range_vectors: + for (unsigned int i = 0; i < m; ++i) + op.reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.reinit_domain_vector = [op](Domain &v, bool fast) + { + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_domain_vectors: + for (unsigned int i = 0; i < m; ++i) + op.reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + return_op.vmult = [op](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + op.vmult(v.block(i), u.block(i)); + }; + + return_op.vmult_add = [op](Range &v, const Domain &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + op.vmult_add(v.block(i), u.block(i)); + }; + + return_op.Tvmult = [op](Domain &v, const Range &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + op.Tvmult(v.block(i), u.block(i)); + }; + + return_op.Tvmult_add = [op](Domain &v, const Range &u) + { + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (unsigned int i = 0; i < m; ++i) + op.Tvmult_add(v.block(i), u.block(i)); + }; + + return return_op; +} + +//@} DEAL_II_NAMESPACE_CLOSE -- 2.39.5