*/
std::function<void(Domain &v, bool fast)> reinit_domain_vector;
+ /**
+ * @name In-place vector space operations
+ */
+ //@{
/**
* Addition with a LinearOperator @p second_op with the same @p Domain
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)$
/**
- * \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)$
/**
- * \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 <typename Range, typename Intermediate, typename Domain>
-LinearOperator<Range, Domain>
-operator*(const LinearOperator<Range, Intermediate> &first_op,
- const LinearOperator<Intermediate, Domain> &second_op)
-{
- LinearOperator<Range, Domain> 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<Intermediate> 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<Intermediate> 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<Intermediate> 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<Intermediate> 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.
/**
- * \relates LinearOperator
+ * @relates LinearOperator
*
* Scalar multiplication of a ScalarOperator object from the right.
*
return number * op;
}
+//@}
+
/**
- * \relates LinearOperator
- *
- * Returns the transpose linear operations of @ref op.
- *
- * @ingroup LAOperators
+ * @name Composition and manipulation of a LinearOperator
*/
-template <typename Range, typename Domain>
-LinearOperator<Domain, Range>
-transpose_operator(const LinearOperator<Range, Domain> &op)
-{
- LinearOperator<Domain, Range> 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 <code>std::function</code> object @ref
- * reinit_vector as an argument to initialize the
- * <code>reinit_range_vector</code> and <code>reinit_domain_vector</code>
- * 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 <typename Range>
-LinearOperator<Range, Range>
-identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
+template <typename Range, typename Intermediate, typename Domain>
+LinearOperator<Range, Domain>
+operator*(const LinearOperator<Range, Intermediate> &first_op,
+ const LinearOperator<Intermediate, Domain> &second_op)
{
- LinearOperator<Range, Range> return_op;
+ LinearOperator<Range, Domain> 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<Intermediate> 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<Intermediate> 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<Intermediate> 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<Intermediate> 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 <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+transpose_operator(const LinearOperator<Range, Domain> &op)
+{
+ LinearOperator<Domain, Range> 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.
*
return return_op;
}
+//@}
+
+
+/**
+ * @name Creation of a LinearOperator
+ */
+//@{
/**
- * \relates LinearOperator
+ * @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.
+ * Returns a LinearOperator that is the identity of the vector space
+ * @p Range.
*
- * The list @p ops is best passed as an initializer list. Consider for
- * example a linear operator block (acting on Vector<double>)
- * @code
- * op_a00 | op_a01
- * |
- * ---------------
- * |
- * op_a10 | op_a11
- * @endcode
- * The coresponding block_operator invocation takes the form
- * @code
- * block_operator<2, 2, BlockVector<double>>({op_a00, op_a01, op_a10, op_a11});
- * @endcode
+ * The function takes an <code>std::function</code> object @ref
+ * reinit_vector as an argument to initialize the
+ * <code>reinit_range_vector</code> and <code>reinit_domain_vector</code>
+ * objects of the LinearOperator object.
*
* @ingroup LAOperators
*/
-template <unsigned int m, unsigned int n,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
-block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+template <typename Range>
+LinearOperator<Range, Range>
+identity_operator(const std::function<void(Range &, bool)> &reinit_vector)
{
- static_assert(m > 0 && n > 0,
- "a blocked LinearOperator must consist of at least one block");
+ LinearOperator<Range, Range> return_op;
- LinearOperator<Range, Domain> 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<double>)
- * <code>diag(op_a0, op_a1, ..., op_am)</code>. The coresponding
- * block_operator invocation takes the form
- * @code
- * block_diagonal_operator<m, BlockVector<double>>({op_00, op_a1, ..., op_am});
- * @endcode
- *
- * @ingroup LAOperators
- */
-template <unsigned int m,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
-block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
-{
- static_assert(m > 0,
- "a blockdiagonal LinearOperator must consist of at least one block");
-
- LinearOperator<Range, Domain> 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 <unsigned int m,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
-block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
-{
-
- static_assert(m > 0,
- "a blockdiagonal LinearOperator must consist of at least "
- "one block");
-
- LinearOperator<Range, Domain> 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
/**
- * \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
/**
- * \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
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<double>)
+ * @code
+ * op_a00 | op_a01
+ * |
+ * ---------------
+ * |
+ * op_a10 | op_a11
+ * @endcode
+ * The coresponding block_operator invocation takes the form
+ * @code
+ * block_operator<2, 2, BlockVector<double>>({op_a00, op_a01, op_a10, op_a11});
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
+template <unsigned int m, unsigned int n,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+LinearOperator<Range, Domain>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &ops)
+{
+ static_assert(m > 0 && n > 0,
+ "a blocked LinearOperator must consist of at least one block");
+
+ LinearOperator<Range, Domain> 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<double>)
+ * <code>diag(op_a0, op_a1, ..., op_am)</code>. The coresponding
+ * block_operator invocation takes the form
+ * @code
+ * block_diagonal_operator<m, BlockVector<double>>({op_00, op_a1, ..., op_am});
+ * @endcode
+ *
+ * @ingroup LAOperators
+ */
+template <unsigned int m,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+LinearOperator<Range, Domain>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &ops)
+{
+ static_assert(m > 0,
+ "a blockdiagonal LinearOperator must consist of at least one block");
+
+ LinearOperator<Range, Domain> 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 <unsigned int m,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+LinearOperator<Range, Domain>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &op)
+{
+
+ static_assert(m > 0,
+ "a blockdiagonal LinearOperator must consist of at least "
+ "one block");
+
+ LinearOperator<Range, Domain> 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
projects / dealii.git / commitdiff