DEAL_II_NAMESPACE_OPEN
// Forward declarations:
+
template <typename Number> class BlockVector;
+template <typename Range = BlockVector<double>,
+ typename Domain = Range>
+class BlockLinearOperator;
+
+template <typename Range = BlockVector<double>,
+ typename Domain = Range,
+ typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_operator(const BlockMatrix &matrix);
+
+template <size_t m, size_t n,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_operator(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, n>, m> &);
+
+template <size_t m,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const std::array<LinearOperator<typename Range::BlockType, typename Domain::BlockType>, m> &);
+
+template <size_t m,
+ typename Range = BlockVector<double>,
+ typename Domain = Range>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const LinearOperator<typename Range::BlockType, typename Domain::BlockType> &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 Range = BlockVector<double>,
+ typename Domain = Range,
+ typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const BlockMatrix &block_matrix);
+
+template <typename Range = BlockVector<double>,
+ typename Domain = Range>
+LinearOperator<Domain, Range>
+block_forward_substitution(const BlockLinearOperator<Range, Domain> &,
+ const BlockLinearOperator<Domain, Range> &);
+
+template <typename Range = BlockVector<double>,
+ typename Domain = Range>
+LinearOperator<Domain, Range>
+block_back_substitution(const BlockLinearOperator<Range, Domain> &,
+ const BlockLinearOperator<Domain, Range> &);
+
+// end of workaround
+
+
+
+/**
+ * TODO: Description
+ */
+template <typename Range, typename Domain>
+class BlockLinearOperator : public LinearOperator<Range, Domain>
+{
+public:
+
+ typedef LinearOperator<typename Range::BlockType, typename Domain::BlockType> BlockType;
+
+ /**
+ * Create an empty BlockLinearOperator object.
+ *
+ * All<code>std::function</code> member objects of this class and its
+ * base class LinearOperator are initialized with default variants that
+ * throw an exception upon invocation.
+ */
+ BlockLinearOperator()
+ : LinearOperator<Range, Domain>()
+ {
+
+ n_block_rows = []() -> unsigned int
+ {
+ Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
+ return 0;
+ };
+
+ n_block_cols = []() -> unsigned int
+ {
+ Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
+ return 0;
+ };
+
+ block = [](unsigned int, unsigned int) -> BlockType
+ {
+ Assert(false, ExcMessage("Uninitialized BlockLinearOperator<Range, Domain>::block called"));
+ return BlockType();
+ };
+ }
+
+ /**
+ * Default copy constructor.
+ */
+ BlockLinearOperator(const BlockLinearOperator<Range, Domain> &) =
+ default;
+
+ /**
+ * Templated copy constructor that creates a BlockLinearOperator object
+ * from an object @p op for which the conversion function
+ * <code>block_operator</code> is defined.
+ */
+ template<typename Op>
+ BlockLinearOperator(const Op &op)
+ {
+ *this = block_operator<Range, Domain, Op>(op);
+ }
+
+ /**
+ * Create a BlockLinearOperator from a two-dimensional array @p ops of
+ * LinearOperator. This constructor calls the corresponding
+ * block_operator() specialization.
+ */
+ template<size_t m, size_t n>
+ BlockLinearOperator(const std::array<std::array<BlockType, n>, m> &ops)
+ {
+ *this = block_operator<m, n, Range, Domain>(ops);
+ }
+
+ /**
+ * Create a block-diagonal BlockLinearOperator from a one-dimensional
+ * array @p ops of LinearOperator. This constructor calls the
+ * corresponding block_operator() specialization.
+ */
+ template<size_t m>
+ BlockLinearOperator(const std::array<BlockType, m> &ops)
+ {
+ *this = block_diagonal_operator<m, Range, Domain>(ops);
+ }
+
+ /**
+ * Default copy assignment operator.
+ */
+ BlockLinearOperator<Range, Domain> &
+ operator=(const BlockLinearOperator<Range, Domain> &) = default;
+
+ /**
+ * Templated copy assignment operator for an object @p op for which the
+ * conversion function <code>block_operator</code> is defined.
+ */
+ template <typename Op>
+ BlockLinearOperator<Range, Domain> &operator=(const Op &op)
+ {
+ *this = block_operator<Range, Domain, Op>(op);
+ return *this;
+ }
+
+ /**
+ * Copy assignment from a two-dimensional array @p ops of LinearOperator.
+ * This assignment operator calls the corresponding block_operator()
+ * specialization.
+ */
+ template <size_t m, size_t n>
+ BlockLinearOperator<Range, Domain> &
+ operator=(const std::array<std::array<BlockType, n>, m> &ops)
+ {
+ *this = block_operator<m, n, Range, Domain>(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 <size_t m>
+ BlockLinearOperator<Range, Domain> &
+ operator=(const std::array<BlockType, m> &ops)
+ {
+ *this = block_diagonal_operator<m, Range, Domain>(ops);
+ return *this;
+ }
+
+ /**
+ * TODO: Description
+ */
+ std::function<unsigned int()> n_block_rows;
+
+ /**
+ * TODO: Description
+ */
+ std::function<unsigned int()> n_block_cols;
+
+ /**
+ * TODO: Description
+ */
+ std::function<BlockType(unsigned int, unsigned int)> block;
+
+ //@}
+};
+
+
+
+namespace internal
+{
+ namespace BlockLinearOperator
+ {
+ // Populate the LinearOperator interfaces with the help of the
+ // BlockLinearOperator functions
+ template <typename Range, typename Domain>
+ inline void
+ populate_linear_operator_functions(
+ dealii::BlockLinearOperator<Range, Domain> &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 <typename Range,
+ typename Domain,
+ typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_operator(const BlockMatrix &block_matrix)
+{
+ typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+
+ BlockLinearOperator<Range, Domain> 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<double>)
*
* @ingroup LAOperators
*/
-template <unsigned int m, unsigned int n,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, size_t n, typename Range, typename Domain>
+BlockLinearOperator<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);
+ typedef typename BlockLinearOperator<Range, Domain>::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<Range, Domain> 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 <typename Range,
+ typename Domain,
+ typename BlockMatrix>
+BlockLinearOperator<Range, Domain>
+block_diagonal_operator(const BlockMatrix &block_matrix)
+{
+ typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+
+ BlockLinearOperator<Range, Domain> 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
*
*
* @ingroup LAOperators
*/
-template <unsigned int m,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, typename Range, typename Domain>
+BlockLinearOperator<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;
+ typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+
+ std::array<std::array<BlockType, m>, 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
*
*
* @ingroup LAOperators
*/
-template <unsigned int m,
- typename Range = BlockVector<double>,
- typename Domain = Range>
-LinearOperator<Range, Domain>
+template <size_t m, typename Range, typename Domain>
+BlockLinearOperator<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));
+ typedef typename BlockLinearOperator<Range, Domain>::BlockType BlockType;
+ std::array<BlockType, m> 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:
*
* and therefore:
* xn = Ann^-1 ( yn - A0n x0 - ... - A(n-1)n x(n-1) )
*
- * @note Notice that we are not using the whole matrix:
- * just the lower triangular block matrix obtained
- * from @p block_matrix it is used.
+ * @note We are not using all blocks of the BlockLinearOperator arguments:
+ * Just the lower triangular block matrix of @p block_operator is used as
+ * well as the diagonal of @p diagonal_inverse.
*
* @ingroup LAOperators
*/
-template < unsigned int n,
- typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
-
-template < unsigned int n,
- typename Range>
-LinearOperator<Range, Range>
-block_forward_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+block_forward_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
+ const BlockLinearOperator<Domain, Range> &diagonal_inverse)
{
LinearOperator<Range, Range> 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<typename Range::BlockType> 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<n; ++i)
+ diagonal_inverse.block(0, 0).vmult(v.block(0), u.block(0));
+ for (unsigned int i = 1; i < m; ++i)
{
- inverse_diagonal[i].reinit_range_vector(*tmp, /*bool fast=*/ true);
- *tmp = u.block(i);
- *tmp *= -1.;
- for (unsigned int j=0; j<i; ++j)
- block_matrix[i][j].vmult_add(*tmp, v.block(j));
- *tmp *= -1.;
- inverse_diagonal[i].vmult(v.block(i),*tmp);
+ auto &dst = v.block(i);
+ dst = u.block(i);
+ dst *= -1.;
+ for (unsigned int j = 0; j < i; ++j)
+ block_operator.block(i, j).vmult_add(dst, v.block(j));
+ dst *= -1.;
+ diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
}
-
- vector_memory.free(tmp);
};
- return_op.vmult_add = [block_matrix, inverse_diagonal](Range &v, const Range &u)
+ return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
{
+ 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));
+
+ if (m == 0)
+ return;
+
static GrowingVectorMemory<typename Range::BlockType> 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<n; ++i)
+ for (unsigned int i = 1; i < m; ++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 (unsigned int j=0; j<i; ++j)
- block_matrix[i][j].vmult_add(*tmp, v.block(j));
+ for (unsigned int j = 0; j < i; ++j)
+ block_operator.block(i, j).vmult_add(*tmp, v.block(j));
*tmp *= -1.;
- inverse_diagonal[i].vmult_add(v.block(i),*tmp);
+ diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
}
vector_memory.free(tmp);
/**
* @relates LinearOperator
*
- * This function implements a back substitution argument to invert an upper
- * block triangular matrix.
- * As argument, it takes an array of array of LinearOperators @p block_matrix
- * representing an upper block triangular matrix and an array of LinearOperators
- * @p inverse_diagonal representing inverses of diagonal blocks of @p block_matrix.
+ * This function implements back substitution to invert an upper block
+ * triangular matrix. As arguments, it takes a BlockLinearOperator
+ * @p block_operator representing an upper block 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:
*
* and therefore:
* x0 = A00^-1 ( y0 - A0n xn - ... - A01 x1 )
*
- * @note Notice that we are not using the whole matrix:
- * just the upper triangular block matrix obtained
- * from @p block_matrix it is used.
+ * @note We are not using all blocks of the BlockLinearOperator arguments:
+ * Just the upper triangular block matrix of @p block_operator is used as
+ * well as the diagonal of @p diagonal_inverse.
+
*
* @ingroup LAOperators
*/
-template < unsigned int n,
- typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &);
-
-template < unsigned int n,
- typename Range>
-LinearOperator<Range, Range>
-block_back_substitution(const std::array<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> &block_matrix,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &inverse_diagonal)
+template <typename Range, typename Domain>
+LinearOperator<Domain, Range>
+block_back_substitution(const BlockLinearOperator<Range, Domain> &block_operator,
+ const BlockLinearOperator<Domain, Range> &diagonal_inverse)
{
LinearOperator<Range, Range> 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<typename Range::BlockType> 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<n; ++j)
- block_matrix[i][j].vmult_add(*tmp,v.block(j));
- *tmp *= -1.;
- inverse_diagonal[i].vmult(v.block(i),*tmp);
+ auto &dst = v.block(i);
+ dst = u.block(i);
+ dst *= -1.;
+ for (int j = i + 1; j < m; ++j)
+ block_operator.block(i, j).vmult_add(dst, v.block(j));
+ dst *= -1.;
+ diagonal_inverse.block(i, i).vmult(dst, dst); // uses intermediate storage
}
-
- vector_memory.free(tmp);
};
- return_op.vmult_add = [&block_matrix, inverse_diagonal](Range &v, const Range &u)
+ return_op.vmult_add = [block_operator, diagonal_inverse](Range &v, const Range &u)
{
+ 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));
static GrowingVectorMemory<typename Range::BlockType> 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<n; ++j)
- block_matrix[i][j].vmult_add(*tmp,v.block(j));
+ for (int j = i + 1; j < m; ++j)
+ block_operator.block(i, j).vmult_add(*tmp,v.block(j));
*tmp *= -1.;
- inverse_diagonal[i].vmult_add(v.block(i),*tmp);
+ diagonal_inverse.block(i, i).vmult_add(v.block(i),*tmp);
}
vector_memory.free(tmp);
return return_op;
}
-/**
- * @relates LinearOperator
- *
- * This function uses above functions block_back_substitution and
- * block_forward_substitution to invert triangular matrices. It takes as input
- * a triangular block matrix @p block_matrix, an array of LinearOperators @p
- * inverse_diagonal representing inverses of block_matrix, and an optional
- * bool @p lower used to specify if block_matrix should be considered as lower
- * triangular matrix (true) or as upper triangular matrix (false). @p lower is
- * equal to true by default.
- *
- */
-
-// 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
-
-template <unsigned int n,
- typename BlockMatrix,
- typename Range = BlockVector<double>>
-LinearOperator<Range, Range>
-block_triangular_inverse(const BlockMatrix &,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n> &,
- bool lower = true);
-
-template <unsigned int n,
- typename BlockMatrix,
- typename Range>
-LinearOperator<Range, Range>
-block_triangular_inverse(const BlockMatrix &block_matrix,
- const std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, 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<std::array<LinearOperator<typename Range::BlockType, typename Range::BlockType>, n>, n> M;
- for (unsigned int i = 0; i<n; ++i)
- for (unsigned int j = 0; j<n; ++j)
- M[i][j] = linear_operator<typename Range::BlockType, typename Range::BlockType>(block_matrix.block(i,j));
-
- if (lower)
- return block_forward_substitution<n, Range>(M, inverse_diagonal);
- else
- return block_back_substitution<n, Range>(M, inverse_diagonal);
-}
-
//@}
DEAL_II_NAMESPACE_CLOSE
//
// ---------------------------------------------------------------------
-// Test block_back_substitution and block_forward_substitution:
+// Test block_operator's constructors and copy assignment variants
#include "../tests.h"
-#include <deal.II/lac/block_linear_operator.h>
-#include <deal.II/lac/block_sparse_matrix.h>
-#include <deal.II/lac/block_vector.h>
-#include <deal.II/lac/dynamic_sparsity_pattern.h>
-#include <deal.II/lac/sparse_matrix.h>
#include <deal.II/lac/vector.h>
-
-#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 <deal.II/lac/block_vector.h>
+#include <deal.II/lac/block_sparse_matrix.h>
+#include <deal.II/lac/block_linear_operator.h>
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<double> 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<double> 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<double> >(a);
-
- auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
- auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
- auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
- auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
- auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
- auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
- auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
- auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
- auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
-
- auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
- auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
- auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
-
- auto inverse_op_a = block_forward_substitution< 3, BlockVector<double> >(
- {{
- {{a00, a01, a02}},
- {{a10, a11, a12}},
- {{a20, a21, a22}}
- }},
- { {d00, d11, d22}});
-
- auto identity = inverse_op_a * op_a;
-
- BlockVector<double> u;
- BlockVector<double> 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<double> a;
+
+ BlockLinearOperator<> op_b(a);
+ auto op_c = block_operator<BlockVector<double>>(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<double> 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<double> 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<double> >(a);
-
- auto a00 = linear_operator< Vector<double>, Vector<double> >(a.block(0,0));
- auto a01 = linear_operator< Vector<double>, Vector<double> >(a.block(0,1));
- auto a02 = linear_operator< Vector<double>, Vector<double> >(a.block(0,2));
- auto a10 = linear_operator< Vector<double>, Vector<double> >(a.block(1,0));
- auto a11 = linear_operator< Vector<double>, Vector<double> >(a.block(1,1));
- auto a12 = linear_operator< Vector<double>, Vector<double> >(a.block(1,2));
- auto a20 = linear_operator< Vector<double>, Vector<double> >(a.block(2,0));
- auto a21 = linear_operator< Vector<double>, Vector<double> >(a.block(2,1));
- auto a22 = linear_operator< Vector<double>, Vector<double> >(a.block(2,2));
-
- auto d00 = linear_operator< Vector<double>, Vector<double> >(d.block(0,0));
- auto d11 = linear_operator< Vector<double>, Vector<double> >(d.block(1,1));
- auto d22 = linear_operator< Vector<double>, Vector<double> >(d.block(2,2));
-
- auto inverse_op_a = block_back_substitution< 3, BlockVector<double> >(
- {{
- {{a00, a01, a02}},
- {{a10, a11, a12}},
- {{a20, a21, a22}}
- }},
- { {d00, d11, d22}});
-
- auto identity = inverse_op_a * op_a;
-
- BlockVector<double> u;
- BlockVector<double> 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<std::array<decltype(op_a), 2>, 2> temp {op_a, op_a, op_a, op_a};
+
+ BlockLinearOperator<> op_b(temp);
+ auto op_c = block_operator<2, 2, BlockVector<double>>(temp);
+ op_c = temp;
+
+ auto op_d = block_operator(temp);
+ op_d = temp;
}
+
+ {
+ std::array<decltype(op_a), 2> temp {op_a, op_a};
+
+ auto op_c = block_diagonal_operator<2, BlockVector<double>>(temp);
+ op_c = temp;
+
+ auto op_d = block_diagonal_operator(temp);
+ op_d = temp;
+ }
+
+ deallog << "OK!" << std::endl;
}
+
projects / dealii.git / commitdiff