From b8fa7b12ef6ec018a11320b13dc20fbf776caaf2 Mon Sep 17 00:00:00 2001 From: Matthias Maier Date: Fri, 7 Aug 2015 21:19:08 -0500 Subject: [PATCH] A BlockLinearOperator class for storing linear operator block structures --- include/deal.II/lac/block_linear_operator.h | 937 +++++++++++------- tests/lac/block_linear_operator_02.cc | 135 +-- tests/lac/block_linear_operator_03.cc | 250 +---- ...ck_linear_operator_03.with_cxx11=on.output | 79 +- tests/lac/block_linear_operator_04.cc | 26 +- 5 files changed, 684 insertions(+), 743 deletions(-) diff --git a/include/deal.II/lac/block_linear_operator.h b/include/deal.II/lac/block_linear_operator.h index 61c84c22bd..690b132a0e 100644 --- a/include/deal.II/lac/block_linear_operator.h +++ b/include/deal.II/lac/block_linear_operator.h @@ -26,23 +26,371 @@ DEAL_II_NAMESPACE_OPEN // Forward declarations: + template class BlockVector; +template , + typename Domain = Range> +class BlockLinearOperator; + +template , + typename Domain = Range, + typename BlockMatrix> +BlockLinearOperator +block_operator(const BlockMatrix &matrix); + +template , + typename Domain = Range> +BlockLinearOperator +block_operator(const std::array, n>, m> &); + +template , + typename Domain = Range> +BlockLinearOperator +block_diagonal_operator(const std::array, m> &); + +template , + typename Domain = Range> +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] +// +// [1] https://gcc.gnu.org/bugzilla/show_bug.cgi?id=53624 +// +// Forward declare functions with partial template defaults: + +template , + typename Domain = Range, + typename BlockMatrix> +BlockLinearOperator +block_diagonal_operator(const BlockMatrix &block_matrix); + +template , + typename Domain = Range> +LinearOperator +block_forward_substitution(const BlockLinearOperator &, + const BlockLinearOperator &); + +template , + typename Domain = Range> +LinearOperator +block_back_substitution(const BlockLinearOperator &, + const BlockLinearOperator &); + +// end of workaround + + + +/** + * TODO: Description + */ +template +class BlockLinearOperator : public LinearOperator +{ +public: + + typedef LinearOperator BlockType; + + /** + * Create an empty BlockLinearOperator object. + * + * Allstd::function member objects of this class and its + * base class LinearOperator are initialized with default variants that + * throw an exception upon invocation. + */ + BlockLinearOperator() + : LinearOperator() + { + + n_block_rows = []() -> unsigned int + { + Assert(false, ExcMessage("Uninitialized BlockLinearOperator::n_block_rows called")); + return 0; + }; + + n_block_cols = []() -> unsigned int + { + Assert(false, ExcMessage("Uninitialized BlockLinearOperator::n_block_cols called")); + return 0; + }; + + block = [](unsigned int, unsigned int) -> BlockType + { + Assert(false, ExcMessage("Uninitialized BlockLinearOperator::block called")); + return BlockType(); + }; + } + + /** + * Default copy constructor. + */ + BlockLinearOperator(const BlockLinearOperator &) = + default; + + /** + * Templated copy constructor that creates a BlockLinearOperator object + * from an object @p op for which the conversion function + * block_operator is defined. + */ + template + BlockLinearOperator(const Op &op) + { + *this = block_operator(op); + } + + /** + * Create a BlockLinearOperator from a two-dimensional array @p ops of + * LinearOperator. This constructor calls the corresponding + * block_operator() specialization. + */ + template + BlockLinearOperator(const std::array, m> &ops) + { + *this = block_operator(ops); + } + + /** + * Create a block-diagonal BlockLinearOperator from a one-dimensional + * array @p ops of LinearOperator. This constructor calls the + * corresponding block_operator() specialization. + */ + template + BlockLinearOperator(const std::array &ops) + { + *this = block_diagonal_operator(ops); + } + + /** + * Default copy assignment operator. + */ + 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) + { + *this = block_operator(op); + return *this; + } + + /** + * Copy assignment from a two-dimensional array @p ops of LinearOperator. + * This assignment operator calls the corresponding block_operator() + * specialization. + */ + template + BlockLinearOperator & + operator=(const std::array, m> &ops) + { + *this = block_operator(ops); + return *this; + } + + /** + * Copy assignment from a one-dimensional array @p ops of LinearOperator + * that creates a block-diagonal BlockLinearOperator. + * This assignment operator calls the corresponding block_operator() + * specialization. + */ + template + BlockLinearOperator & + operator=(const std::array &ops) + { + *this = block_diagonal_operator(ops); + return *this; + } + + /** + * TODO: Description + */ + std::function n_block_rows; + + /** + * TODO: Description + */ + std::function n_block_cols; + + /** + * TODO: Description + */ + std::function block; + + //@} +}; + + + +namespace internal +{ + namespace BlockLinearOperator + { + // Populate the LinearOperator interfaces with the help of the + // BlockLinearOperator functions + template + inline void + populate_linear_operator_functions( + dealii::BlockLinearOperator &op) + { + op.reinit_range_vector = [=](Range &v, bool fast) + { + const unsigned int m = op.n_block_rows(); + + // Reinitialize the block vector to m blocks: + v.reinit(m); + + // And reinitialize every individual block with reinit_range_vectors: + for (size_t i = 0; i < m; ++i) + op.block(i, 0).reinit_range_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + op.reinit_domain_vector = [=](Domain &v, bool fast) + { + const unsigned int n = op.n_block_cols(); + + // Reinitialize the block vector to n blocks: + v.reinit(n); + + // And reinitialize every individual block with reinit_domain_vectors: + for (size_t i = 0; i < n; ++i) + op.block(0, i).reinit_domain_vector(v.block(i), fast); + + v.collect_sizes(); + }; + + op.vmult = [=](Range &v, const Domain &u) + { + const unsigned int m = op.n_block_rows(); + const unsigned int n = op.n_block_cols(); + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); + + for (size_t i = 0; i < m; ++i) + { + op.block(i, 0).vmult(v.block(i), u.block(0)); + for (size_t j = 1; j < n; ++j) + op.block(i, j).vmult_add(v.block(i), u.block(j)); + } + }; + + op.vmult_add = [=](Range &v, const Domain &u) + { + const unsigned int m = op.n_block_rows(); + const unsigned int n = op.n_block_cols(); + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); + + for (size_t i = 0; i < m; ++i) + for (size_t j = 0; j < n; ++j) + op.block(i, j).vmult_add(v.block(i), u.block(j)); + }; + + op.Tvmult = [=](Domain &v, const Range &u) + { + const unsigned int n = op.n_block_cols(); + const unsigned int m = op.n_block_rows(); + Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (size_t i = 0; i < n; ++i) + { + op.block(0, i).Tvmult(v.block(i), u.block(0)); + for (size_t j = 1; j < m; ++j) + op.block(j, i).Tvmult_add(v.block(i), u.block(j)); + } + }; + + op.Tvmult_add = [=](Domain &v, const Range &u) + { + const unsigned int n = op.n_block_cols(); + const unsigned int m = op.n_block_rows(); + Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + + for (size_t i = 0; i < n; ++i) + for (size_t j = 0; j < m; ++j) + op.block(j, i).Tvmult_add(v.block(i), u.block(j)); + }; + } + } /*namespace BlockLinearOperator*/ +} /*namespace internal*/ + + + /** - * @name Creation of LinearOperator block structures + * @name Creation of a BlockLinearOperator */ //@{ /** * @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. + * A function that encapsulates a @p block_matrix into a + * BlockLinearOperator. + * + * All changes made on the block structure and individual blocks of @p + * block_matrix after the creation of the BlockLinearOperator object are + * reflected by the operator object. + * + * @ingroup LAOperators + */ +template +BlockLinearOperator +block_operator(const BlockMatrix &block_matrix) +{ + typedef typename BlockLinearOperator::BlockType BlockType; + + BlockLinearOperator return_op; + + return_op.n_block_rows = [&block_matrix]() -> unsigned int + { + return block_matrix.n_block_rows(); + }; + + return_op.n_block_cols = [&block_matrix]() -> unsigned int + { + return block_matrix.n_block_cols(); + }; + + return_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType + { +#ifdef DEBUG + const unsigned int m = block_matrix.n_block_rows(); + const unsigned int n = block_matrix.n_block_cols(); + Assert(i < m, ExcIndexRange (i, 0, m)); + Assert(j < n, ExcIndexRange (j, 0, n)); +#endif + + return BlockType(block_matrix.block(i, j)); + }; + + internal::BlockLinearOperator::populate_linear_operator_functions(return_op); + return return_op; +} + + + +/** + * @relates LinearOperator + * + * A variant of above function that encapsulates a given collection @p ops + * of LinearOperators into a block structure. Here, 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 + * 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) @@ -60,91 +408,97 @@ template class BlockVector; * * @ingroup LAOperators */ -template , - typename Domain = Range> -LinearOperator +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"); - LinearOperator return_op; - - return_op.reinit_range_vector = [ops](Range &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(m); + typedef typename BlockLinearOperator::BlockType BlockType; - // 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); + BlockLinearOperator return_op; - v.collect_sizes(); + return_op.n_block_rows = []() -> unsigned int + { + return m; }; - return_op.reinit_domain_vector = [ops](Domain &v, bool fast) + return_op.n_block_cols = []() -> unsigned int { - // 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 n; }; - return_op.vmult = [ops](Range &v, const Domain &u) + return_op.block = [ops](unsigned int i, unsigned int j) -> BlockType { - Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); - Assert(u.n_blocks() == n, ExcDimensionMismatch(u.n_blocks(), n)); + Assert(i < m, ExcIndexRange (i, 0, m)); + Assert(j < n, ExcIndexRange (j, 0, 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 ops[i][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)); + internal::BlockLinearOperator::populate_linear_operator_functions(return_op); + return return_op; +} - 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)); - } +/** + * @relates LinearOperator + * + * This function extracts the diagonal blocks of @p block_matrix (either a + * block matrix type or a BlockLinearOperator) and creates a + * BlockLinearOperator with the diagonal. Off-diagonal elements are + * initialized as null_operator (with correct reinit_range_vector and + * reinit_domain_vector methods). + * + * All changes made on the individual diagonal blocks of @p block_matrix + * after the creation of the BlockLinearOperator object are reflected by + * the operator object. + * + * @ingroup LAOperators + */ +template +BlockLinearOperator +block_diagonal_operator(const BlockMatrix &block_matrix) +{ + typedef typename BlockLinearOperator::BlockType BlockType; + + BlockLinearOperator return_op; + + return_op.n_block_rows = [&block_matrix]() -> unsigned int + { + return block_matrix.n_block_rows(); }; - return_op.Tvmult_add = [ops](Domain &v, const Range &u) + return_op.n_block_cols = [&block_matrix]() -> unsigned int { - Assert(v.n_blocks() == n, ExcDimensionMismatch(v.n_blocks(), n)); - Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); + return block_matrix.n_block_cols(); + }; - 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_op.block = [&block_matrix](unsigned int i, unsigned int j) -> BlockType + { +#ifdef DEBUG + const unsigned int m = block_matrix.n_block_rows(); + const unsigned int n = block_matrix.n_block_cols(); + Assert(m == n, ExcDimensionMismatch(m, n)); + Assert(i < m, ExcIndexRange (i, 0, m)); + Assert(j < n, ExcIndexRange (j, 0, n)); +#endif + if (i == j) + return BlockType(block_matrix.block(i, j)); + else + return null_operator(BlockType(block_matrix.block(i, j))); }; + internal::BlockLinearOperator::populate_linear_operator_functions(return_op); return return_op; } + /** * @relates LinearOperator * @@ -162,81 +516,41 @@ block_operator(const std::array, - typename Domain = Range> -LinearOperator +template +BlockLinearOperator 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; + typedef typename BlockLinearOperator::BlockType BlockType; + + std::array, m> new_ops; + + // This is a bit tricky. We have to make sure that the off-diagonal + // elements of return_op.ops are populated correctly. They must be + // null_operators, but with correct reinit_domain_vector and + // reinit_range_vector functions. + for (size_t i = 0; i < m; ++i) + for (size_t j = 0; j < m; ++j) + if (i == j) + { + // diagonal elements are easy: + new_ops[i][j] = ops[i]; + } + else + { + // create a null-operator... + new_ops[i][j] = null_operator(ops[i]); + // ... and fix up reinit_domain_vector: + new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector; + } + + return block_operator(new_ops); } + /** * @relates LinearOperator * @@ -246,91 +560,29 @@ block_diagonal_operator(const std::array, - typename Domain = Range> -LinearOperator +template +BlockLinearOperator 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)); + typedef typename BlockLinearOperator::BlockType BlockType; + std::array new_ops; + new_ops.fill(op); - 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; + return block_diagonal_operator(new_ops); } - /** * @relates LinearOperator * - * This function implements a forward substitution argument to invert a lower - * block triangular matrix. - * As argument, it takes an array of array of LinearOperators @p block_matrix - * representing a lower block triangular matrix and an array of LinearOperators - * @p inverse_diagonal representing inverses of diagonal blocks of @p block_matrix. + * This function implements forward substitution to invert a lower block + * triangular matrix. As arguments, it takes a BlockLinearOperator + * @p block_operator representing a block lower triangular matrix, as well as + * a BlockLinearOperator @p diagonal_inverse representing inverses of + * diagonal blocks of @p block_operator. * * Let us assume we have a linear system with the following block structure: * @@ -345,85 +597,79 @@ block_diagonal_operator(const LinearOperator> -LinearOperator -block_forward_substitution(const std::array, n>, n> &, - const std::array, n> &); - -template < unsigned int n, - typename Range> -LinearOperator -block_forward_substitution(const std::array, n>, n> &block_matrix, - const std::array, n> &inverse_diagonal) +template +LinearOperator +block_forward_substitution(const BlockLinearOperator &block_operator, + const BlockLinearOperator &diagonal_inverse) { LinearOperator return_op; - return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(n); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < n; ++i) - inverse_diagonal[i].reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; + return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector; + return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector; - return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast) + return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u) { - // Reinitialize the block vector to m blocks: - v.reinit(n); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < n; ++i) - inverse_diagonal[i].reinit_domain_vector(v.block(i), fast); - - v.collect_sizes(); - }; + const unsigned int m = block_operator.n_block_rows(); + Assert(block_operator.n_block_cols() == m, + ExcDimensionMismatch(block_operator.n_block_cols(), m)); + Assert(diagonal_inverse.n_block_rows() == m, + ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m)); + Assert(diagonal_inverse.n_block_cols() == m, + ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m)); + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u) - { - static GrowingVectorMemory vector_memory; - typename Range::BlockType *tmp = vector_memory.alloc(); + if (m == 0) + return; - inverse_diagonal[0].vmult(v.block(0), u.block(0)); - for (unsigned int i=1; i vector_memory; - typename Range::BlockType *tmp = vector_memory.alloc(); + typename Range::BlockType *tmp = vector_memory.alloc(); + + diagonal_inverse.block(0, 0).vmult_add(v.block(0), u.block(0)); - inverse_diagonal[0].vmult_add(v.block(0), u.block(0)); - for (unsigned int i=1; i> -LinearOperator -block_back_substitution(const std::array, n>, n> &, - const std::array, n> &); - -template < unsigned int n, - typename Range> -LinearOperator -block_back_substitution(const std::array, n>, n> &block_matrix, - const std::array, n> &inverse_diagonal) +template +LinearOperator +block_back_substitution(const BlockLinearOperator &block_operator, + const BlockLinearOperator &diagonal_inverse) { LinearOperator return_op; - return_op.reinit_range_vector = [inverse_diagonal](Range &v, bool fast) - { - // Reinitialize the block vector to m blocks: - v.reinit(n); - - // And reinitialize every individual block with reinit_range_vectors: - for (unsigned int i = 0; i < n; ++i) - inverse_diagonal[i].reinit_range_vector(v.block(i), fast); - - v.collect_sizes(); - }; + return_op.reinit_range_vector = diagonal_inverse.reinit_range_vector; + return_op.reinit_domain_vector = diagonal_inverse.reinit_domain_vector; - return_op.reinit_domain_vector = [inverse_diagonal](Range &v, bool fast) + return_op.vmult = [block_operator, diagonal_inverse](Range &v, const Range &u) { - // Reinitialize the block vector to m blocks: - v.reinit(n); - - // And reinitialize every individual block with reinit_domain_vectors: - for (unsigned int i = 0; i < n; ++i) - inverse_diagonal[i].reinit_domain_vector(v.block(i), fast); + const unsigned int m = block_operator.n_block_rows(); + Assert(block_operator.n_block_cols() == m, + ExcDimensionMismatch(block_operator.n_block_cols(), m)); + Assert(diagonal_inverse.n_block_rows() == m, + ExcDimensionMismatch(diagonal_inverse.n_block_rows(), m)); + Assert(diagonal_inverse.n_block_cols() == m, + ExcDimensionMismatch(diagonal_inverse.n_block_cols(), m)); + Assert(v.n_blocks() == m, ExcDimensionMismatch(v.n_blocks(), m)); + Assert(u.n_blocks() == m, ExcDimensionMismatch(u.n_blocks(), m)); - v.collect_sizes(); - }; + if (m == 0) + return; - return_op.vmult = [block_matrix, inverse_diagonal](Range &v, const Range &u) - { - static GrowingVectorMemory vector_memory; - typename Range::BlockType *tmp = vector_memory.alloc(); + diagonal_inverse.block(m-1, m-1).vmult(v.block(m-1),u.block(m-1)); - inverse_diagonal[n-1].vmult(v.block(n-1),u.block(n-1)); - for (int i=n-2; i>=0; --i) + for (int i = m - 2; i >= 0; --i) { - inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true); - *tmp = u.block(i); - *tmp *= -1.; - for (int j=i+1; j vector_memory; typename Range::BlockType *tmp = vector_memory.alloc(); - inverse_diagonal[n-1].vmult_add(v.block(n-1),u.block(n-1)); - for (int i=n-2; i>=0; --i) + if (m == 0) + return; + + diagonal_inverse.block(m-1, m-1).vmult_add(v.block(m-1),u.block(m-1)); + + for (int i = m - 2; i >= 0; --i) { - inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true); + diagonal_inverse.block(i, i).reinit_range_vector(*tmp, /*bool fast=*/ true); *tmp = u.block(i); *tmp *= -1.; - for (int j=i+1; j> -LinearOperator -block_triangular_inverse(const BlockMatrix &, - const std::array, n> &, - bool lower = true); - -template -LinearOperator -block_triangular_inverse(const BlockMatrix &block_matrix, - const std::array, n> &inverse_diagonal, - bool lower) -{ - Assert(block_matrix.n_block_rows() == n, - ExcDimensionMismatch(block_matrix.n_block_rows(), n)); - Assert(block_matrix.n_block_rows() == block_matrix.n_block_cols(), - ExcDimensionMismatch(block_matrix.n_block_rows(), - block_matrix.n_block_cols())); - - std::array, n>, n> M; - for (unsigned int i = 0; i(block_matrix.block(i,j)); - - if (lower) - return block_forward_substitution(M, inverse_diagonal); - else - return block_back_substitution(M, inverse_diagonal); -} - //@} DEAL_II_NAMESPACE_CLOSE diff --git a/tests/lac/block_linear_operator_02.cc b/tests/lac/block_linear_operator_02.cc index 89b2827819..56b9a614e9 100644 --- a/tests/lac/block_linear_operator_02.cc +++ b/tests/lac/block_linear_operator_02.cc @@ -78,7 +78,7 @@ int main() a.block(2, 2).set(i, i, 3.); - auto op_a = linear_operator>(a); + auto op_a = block_operator(a); auto op_b0 = linear_operator(a.block(0, 0)); auto op_b1 = linear_operator(a.block(1, 1)); @@ -88,68 +88,71 @@ int main() std::array temp{op_b0, op_b1, op_b2}; auto op_b = block_diagonal_operator<3, BlockVector>(temp); - // vmult: - BlockVector u; - op_a.reinit_domain_vector(u, false); - for (unsigned int i = 0; i < u.size(); ++i) { - u[i] = (double)(i+1); - } + { + + BlockVector u; + op_a.reinit_domain_vector(u, false); + for (unsigned int i = 0; i < u.size(); ++i) { + u[i] = (double)(i+1); + } - PRINTME("u", u); + PRINTME("u", u); - BlockVector v; - op_a.reinit_domain_vector(v, false); - BlockVector w; - op_a.reinit_domain_vector(w, false); - BlockVector x; - op_a.reinit_domain_vector(x, false); + BlockVector v; + op_a.reinit_domain_vector(v, false); + BlockVector w; + op_a.reinit_domain_vector(w, false); + BlockVector x; + op_a.reinit_domain_vector(x, false); - op_a.vmult(v, u); - PRINTME("Au", v); + op_a.vmult(v, u); + PRINTME("Au", v); - op_b.vmult(w, u); - PRINTME("Bu", w); + op_b.vmult(w, u); + PRINTME("Bu", w); - x = v; - x -= w; - PRINTME("Au-Bu", x); + x = v; + x -= w; + PRINTME("Au-Bu", x); - // Test that both objects give the same results: + // Test that both objects give the same results: - auto op_x = op_a - op_b; + auto op_x = op_a - op_b; - op_x.vmult(x, u); - PRINTME("(A-B).vmult", x); + op_x.vmult(x, u); + PRINTME("(A-B).vmult", x); - x = 0.; - op_x.vmult_add(x, u); - PRINTME("(A-B).vmult_add", x); + x = 0.; + op_x.vmult_add(x, u); + PRINTME("(A-B).vmult_add", x); - op_x.Tvmult(x, u); - PRINTME("(A-B).Tvmult", x); + op_x.Tvmult(x, u); + PRINTME("(A-B).Tvmult", x); - x = 0.; - op_x.Tvmult_add(x, u); - PRINTME("(A-B).Tvmult_add", x); + x = 0.; + op_x.Tvmult_add(x, u); + PRINTME("(A-B).Tvmult_add", x); - // Test vector reinitalization: + // Test vector reinitalization: - op_x = op_b * op_b * op_b; - op_x.vmult(x, u); - PRINTME("(B*B*B) vmult", x); + op_x = op_b * op_b * op_b; + op_x.vmult(x, u); + PRINTME("(B*B*B) vmult", x); - x = 0.; - op_x.vmult_add(x, u); - PRINTME("(B*B*B) vmult_add", x); + x = 0.; + op_x.vmult_add(x, u); + PRINTME("(B*B*B) vmult_add", x); - op_x.Tvmult(x, u); - PRINTME("(B*B*B) Tvmult", x); + op_x.Tvmult(x, u); + PRINTME("(B*B*B) Tvmult", x); - x = 0.; - op_x.Tvmult_add(x, u); - PRINTME("(B*B*B) Tvmult_add", x); + x = 0.; + op_x.Tvmult_add(x, u); + PRINTME("(B*B*B) Tvmult_add", x); + + } // And finally the other block_diagonal_operator variant: @@ -158,30 +161,32 @@ int main() auto op_d = block_diagonal_operator<5, BlockVector>(op_b0); - op_c.reinit_domain_vector(u, false); - for (unsigned int i = 0; i < u.size(); ++i) { - u[i] = (double)(i+1); - } - PRINTME("u", u); - op_c.reinit_domain_vector(x, false); - - op_x = op_c - op_d; - - op_x.vmult(x, u); - PRINTME("(C-D) vmult", x); + { + BlockVector u; + op_c.reinit_domain_vector(u, false); + for (unsigned int i = 0; i < u.size(); ++i) { + u[i] = (double)(i+1); + } + PRINTME("u", u); - x = 0.; - op_x.vmult_add(x, u); - PRINTME("(C-D) vmult_add", x); + BlockVector x; + op_c.reinit_range_vector(x, false); - op_x.Tvmult(x, u); - PRINTME("(C-D) Tvmult", x); + auto op_x = op_c - op_d; - x = 0.; - op_x.Tvmult_add(x, u); - PRINTME("(C-D) Tvmult_add", x); + op_x.vmult(x, u); + PRINTME("(C-D) vmult", x); -} + x = 0.; + op_x.vmult_add(x, u); + PRINTME("(C-D) vmult_add", x); + op_x.Tvmult(x, u); + PRINTME("(C-D) Tvmult", x); + x = 0.; + op_x.Tvmult_add(x, u); + PRINTME("(C-D) Tvmult_add", x); + } +} diff --git a/tests/lac/block_linear_operator_03.cc b/tests/lac/block_linear_operator_03.cc index 7c8378cb2a..4ec6e5ecf8 100644 --- a/tests/lac/block_linear_operator_03.cc +++ b/tests/lac/block_linear_operator_03.cc @@ -13,237 +13,55 @@ // // --------------------------------------------------------------------- -// Test block_back_substitution and block_forward_substitution: +// Test block_operator's constructors and copy assignment variants #include "../tests.h" -#include -#include -#include -#include -#include #include - -#define PRINTME(name, var) \ - deallog << "Block vector: " name << ":" << std::endl; \ - for (unsigned int i = 0; i < var.n_blocks(); ++i) \ - deallog << "[block " << i << " ] " << var.block(i); - +#include +#include +#include using namespace dealii; int main() { initlog(); - deallog << std::setprecision(12); - // BlockSparseMatrix: { - BlockDynamicSparsityPattern dsp(3, 3); - for (unsigned int i = 0; i < 3; ++i) - for (unsigned int j = 0; j < 3; ++j) - dsp.block(i, j).reinit (1, 1); - dsp.collect_sizes (); - - BlockSparsityPattern sparsity_pattern; - sparsity_pattern.copy_from(dsp); - sparsity_pattern.compress(); - - BlockSparseMatrix a (sparsity_pattern); - for (unsigned int i = 0; i < 3; ++i) - { - a.block(i,i).set(0, 0, i+i +1); - for (unsigned int j = 0; j < i; ++j) - a.block(i,j).set(0, 0, 10); - } - - BlockSparseMatrix d(sparsity_pattern); - for (unsigned int i = 0; i < 3; ++i) - d.block(i,i).set(0,0, 1.0 / (i+i +1) ); - - auto op_a = linear_operator< BlockVector >(a); - - auto a00 = linear_operator< Vector, Vector >(a.block(0,0)); - auto a01 = linear_operator< Vector, Vector >(a.block(0,1)); - auto a02 = linear_operator< Vector, Vector >(a.block(0,2)); - auto a10 = linear_operator< Vector, Vector >(a.block(1,0)); - auto a11 = linear_operator< Vector, Vector >(a.block(1,1)); - auto a12 = linear_operator< Vector, Vector >(a.block(1,2)); - auto a20 = linear_operator< Vector, Vector >(a.block(2,0)); - auto a21 = linear_operator< Vector, Vector >(a.block(2,1)); - auto a22 = linear_operator< Vector, Vector >(a.block(2,2)); - - auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); - auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); - auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); - - auto inverse_op_a = block_forward_substitution< 3, BlockVector >( - {{ - {{a00, a01, a02}}, - {{a10, a11, a12}}, - {{a20, a21, a22}} - }}, - { {d00, d11, d22}}); - - auto identity = inverse_op_a * op_a; - - BlockVector u; - BlockVector v; - - deallog << " -- Matrix -- " << std::endl; - op_a.reinit_domain_vector(u, false); - op_a.reinit_range_vector(v, false); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0;; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1; - - op_a.vmult(v, u); - - PRINTME("v", v); - } - - deallog << " -- Inverse -- " << std::endl; - inverse_op_a.reinit_domain_vector(u, false); - inverse_op_a.reinit_range_vector(v, true); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1;; - - inverse_op_a.vmult(v, u); - - PRINTME("v", v); - } - - deallog << " -- Identity -- " << std::endl; - identity.reinit_domain_vector(u, false); - identity.reinit_range_vector(v, false); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0;; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1;; - - identity.vmult(v, u); - - PRINTME("v", v); - } + BlockSparseMatrix a; + + BlockLinearOperator<> op_b(a); + auto op_c = block_operator>(a); + op_c = a; + + auto op_d = block_operator(a); + op_d = a; } + auto op_a = LinearOperator<>(); { - BlockDynamicSparsityPattern dsp(3, 3); - for (unsigned int i = 0; i < 3; ++i) - for (unsigned int j = 0; j < 3; ++j) - dsp.block(i, j).reinit (1, 1); - dsp.collect_sizes (); - - BlockSparsityPattern sparsity_pattern; - sparsity_pattern.copy_from(dsp); - sparsity_pattern.compress(); - - BlockSparseMatrix a (sparsity_pattern); - for (unsigned int i = 0; i < 3; ++i) - { - a.block(i,i).set(0, 0, i+i +1); - for (unsigned int j = i+1; j < 3; ++j) - a.block(i,j).set(0, 0, 10); - } - - BlockSparseMatrix d(sparsity_pattern); - for (unsigned int i = 0; i < 3; ++i) - d.block(i,i).set(0,0, 1.0 / (i+i +1) ); - - auto op_a = linear_operator< BlockVector >(a); - - auto a00 = linear_operator< Vector, Vector >(a.block(0,0)); - auto a01 = linear_operator< Vector, Vector >(a.block(0,1)); - auto a02 = linear_operator< Vector, Vector >(a.block(0,2)); - auto a10 = linear_operator< Vector, Vector >(a.block(1,0)); - auto a11 = linear_operator< Vector, Vector >(a.block(1,1)); - auto a12 = linear_operator< Vector, Vector >(a.block(1,2)); - auto a20 = linear_operator< Vector, Vector >(a.block(2,0)); - auto a21 = linear_operator< Vector, Vector >(a.block(2,1)); - auto a22 = linear_operator< Vector, Vector >(a.block(2,2)); - - auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); - auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); - auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); - - auto inverse_op_a = block_back_substitution< 3, BlockVector >( - {{ - {{a00, a01, a02}}, - {{a10, a11, a12}}, - {{a20, a21, a22}} - }}, - { {d00, d11, d22}}); - - auto identity = inverse_op_a * op_a; - - BlockVector u; - BlockVector v; - - deallog << " -- Matrix -- " << std::endl; - op_a.reinit_domain_vector(u, false); - op_a.reinit_range_vector(v, false); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0;; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1; - - op_a.vmult(v, u); - - PRINTME("v", v); - } - - deallog << " -- Inverse -- " << std::endl; - inverse_op_a.reinit_domain_vector(u, false); - inverse_op_a.reinit_range_vector(v, true); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1;; - - inverse_op_a.vmult(v, u); - - PRINTME("v", v); - } - - deallog << " -- Identity -- " << std::endl; - identity.reinit_domain_vector(u, false); - identity.reinit_range_vector(v, false); - for(unsigned int j = 0; j<3; ++j) - { - for(unsigned int i = 0; i<3; ++i) - { - u.block(i)[0] = 0;; - v.block(i)[0] = 0; - } - u.block(j)[0] = 1;; - - identity.vmult(v, u); - - PRINTME("v", v); - } + std::array, 2> temp {op_a, op_a, op_a, op_a}; + + BlockLinearOperator<> op_b(temp); + auto op_c = block_operator<2, 2, BlockVector>(temp); + op_c = temp; + + auto op_d = block_operator(temp); + op_d = temp; } + + { + std::array temp {op_a, op_a}; + + auto op_c = block_diagonal_operator<2, BlockVector>(temp); + op_c = temp; + + auto op_d = block_diagonal_operator(temp); + op_d = temp; + } + + deallog << "OK!" << std::endl; } + diff --git a/tests/lac/block_linear_operator_03.with_cxx11=on.output b/tests/lac/block_linear_operator_03.with_cxx11=on.output index 1e80e98cd0..5cfb783b8f 100644 --- a/tests/lac/block_linear_operator_03.with_cxx11=on.output +++ b/tests/lac/block_linear_operator_03.with_cxx11=on.output @@ -1,79 +1,2 @@ -DEAL:: -- Matrix -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] 10.0000000000 -DEAL::[block 2 ] 10.0000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 3.00000000000 -DEAL::[block 2 ] 10.0000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 5.00000000000 -DEAL:: -- Inverse -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] -3.33333333333 -DEAL::[block 2 ] 4.66666666667 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 0.333333333333 -DEAL::[block 2 ] -0.666666666667 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 0.200000000000 -DEAL:: -- Identity -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 1.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 1.00000000000 -DEAL:: -- Matrix -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 10.0000000000 -DEAL::[block 1 ] 3.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 10.0000000000 -DEAL::[block 1 ] 10.0000000000 -DEAL::[block 2 ] 5.00000000000 -DEAL:: -- Inverse -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] -3.33333333333 -DEAL::[block 1 ] 0.333333333333 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 4.66666666667 -DEAL::[block 1 ] -0.666666666667 -DEAL::[block 2 ] 0.200000000000 -DEAL:: -- Identity -- -DEAL::Block vector: v: -DEAL::[block 0 ] 1.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 1.00000000000 -DEAL::[block 2 ] 0.00000000000 -DEAL::Block vector: v: -DEAL::[block 0 ] 0.00000000000 -DEAL::[block 1 ] 0.00000000000 -DEAL::[block 2 ] 1.00000000000 +DEAL::OK! diff --git a/tests/lac/block_linear_operator_04.cc b/tests/lac/block_linear_operator_04.cc index e49548d508..07cbca9e8a 100644 --- a/tests/lac/block_linear_operator_04.cc +++ b/tests/lac/block_linear_operator_04.cc @@ -61,12 +61,14 @@ int main() for (unsigned int i = 0; i < 3; ++i) d.block(i,i).set(0,0, 1.0 / (i+i +1) ); - auto op_a = linear_operator< BlockVector >(a); - auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); - auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); - auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + auto op_a = linear_operator>(a); + + auto op_b1 = block_operator(a); + auto op_b2 = block_diagonal_operator(d); + + auto inverse_op_a = + block_forward_substitution>(op_b1, op_b2); - auto inverse_op_a = block_triangular_inverse<3, BlockSparseMatrix >(a, {{d00, d11, d22}}); auto identity = inverse_op_a * op_a; BlockVector u; @@ -146,14 +148,16 @@ int main() BlockSparseMatrix d(sparsity_pattern); for (unsigned int i = 0; i < 3; ++i) - d.block(i,i).set(0,0, 1.0 / (i+i +1) ); + d.block(i, i).set(0, 0, 1.0 / (i + i + 1)); + + auto op_a = linear_operator>(a); + + auto op_b1 = block_operator(a); + auto op_b2 = block_diagonal_operator(d); - auto op_a = linear_operator< BlockVector >(a); - auto d00 = linear_operator< Vector, Vector >(d.block(0,0)); - auto d11 = linear_operator< Vector, Vector >(d.block(1,1)); - auto d22 = linear_operator< Vector, Vector >(d.block(2,2)); + auto inverse_op_a = + block_back_substitution>(op_b1, op_b2); - auto inverse_op_a = block_triangular_inverse< 3, BlockSparseMatrix >(a, {{d00, d11, d22}}, false); auto identity = inverse_op_a * op_a; BlockVector u; -- 2.39.5