15#ifndef dealii_block_linear_operator_h
16#define dealii_block_linear_operator_h
31 namespace BlockLinearOperatorImplementation
33 template <
typename PayloadBlockType =
35 class EmptyBlockPayload;
39template <
typename Number>
42template <
typename Range = BlockVector<
double>,
43 typename Domain = Range,
44 typename BlockPayload =
45 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
49template <
typename Range = BlockVector<
double>,
50 typename Domain = Range,
51 typename BlockPayload =
52 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
53 typename BlockMatrixType>
57template <std::size_t m,
60 typename Domain = Range,
61 typename BlockPayload =
65 const std::array<std::array<
LinearOperator<
typename Range::BlockType,
66 typename Domain::BlockType,
67 typename BlockPayload::BlockType>,
71template <std::size_t m,
73 typename Domain = Range,
74 typename BlockPayload =
79 typename Domain::BlockType,
80 typename BlockPayload::BlockType>,
83template <std::size_t m,
85 typename Domain = Range,
86 typename BlockPayload =
91 typename Domain::BlockType,
92 typename BlockPayload::BlockType> &op);
164template <
typename Range,
typename Domain,
typename BlockPayload>
166 :
public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
170 typename Domain::BlockType,
171 typename BlockPayload::BlockType>;
182 typename BlockPayload::
BlockType(payload, payload))
188 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
196 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
204 "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
220 template <
typename Op>
223 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
231 template <std::
size_t m, std::
size_t n>
234 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
242 template <std::
size_t m>
245 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
258 template <
typename Op>
262 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
271 template <std::
size_t m, std::
size_t n>
273 operator=(
const std::array<std::array<BlockType, n>, m> &ops)
275 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
273 operator=(
const std::array<std::array<BlockType, n>, m> &ops) {
…}
284 template <std::
size_t m>
288 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
315 namespace BlockLinearOperatorImplementation
322 template <
typename Function1,
328 const Function2 &loop_op,
336 tmp->reinit(v,
true);
338 const unsigned int n = u.n_blocks();
339 const unsigned int m = v.n_blocks();
341 for (
unsigned int i = 0; i < m; ++i)
343 first_op(*tmp, u, i, 0);
344 for (
unsigned int j = 1; j < n; ++j)
345 loop_op(*tmp, u, i, j);
356 template <
typename Range,
typename Domain,
typename BlockPayload>
368 for (
unsigned int i = 0; i < m; ++i)
369 op.
block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
381 for (
unsigned int i = 0; i < n; ++i)
382 op.
block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
387 op.
vmult = [&op](Range &v,
const Domain &u) {
395 const auto first_op = [&op](Range &v,
397 const unsigned int i,
398 const unsigned int j) {
399 op.
block(i, j).vmult(v.block(i), u.block(j));
402 const auto loop_op = [&op](Range &v,
404 const unsigned int i,
405 const unsigned int j) {
406 op.
block(i, j).vmult_add(v.block(i), u.block(j));
413 for (
unsigned int i = 0; i < m; ++i)
415 op.
block(i, 0).vmult(v.block(i), u.block(0));
416 for (
unsigned int j = 1; j < n; ++j)
417 op.
block(i, j).vmult_add(v.block(i), u.block(j));
422 op.
vmult_add = [&op](Range &v,
const Domain &u) {
430 const auto first_op = [&op](Range &v,
432 const unsigned int i,
433 const unsigned int j) {
434 op.
block(i, j).vmult(v.block(i), u.block(j));
437 const auto loop_op = [&op](Range &v,
439 const unsigned int i,
440 const unsigned int j) {
441 op.
block(i, j).vmult_add(v.block(i), u.block(j));
448 for (
unsigned int i = 0; i < m; ++i)
449 for (
unsigned int j = 0; j < n; ++j)
450 op.
block(i, j).vmult_add(v.block(i), u.block(j));
454 op.
Tvmult = [&op](Domain &v,
const Range &u) {
462 const auto first_op = [&op](Range &v,
464 const unsigned int i,
465 const unsigned int j) {
466 op.
block(j, i).Tvmult(v.block(i), u.block(j));
469 const auto loop_op = [&op](Range &v,
471 const unsigned int i,
472 const unsigned int j) {
473 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
480 for (
unsigned int i = 0; i < n; ++i)
482 op.
block(0, i).Tvmult(v.block(i), u.block(0));
483 for (
unsigned int j = 1; j < m; ++j)
484 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
489 op.
Tvmult_add = [&op](Domain &v,
const Range &u) {
497 const auto first_op = [&op](Range &v,
499 const unsigned int i,
500 const unsigned int j) {
501 op.
block(j, i).Tvmult(v.block(i), u.block(j));
504 const auto loop_op = [&op](Range &v,
506 const unsigned int i,
507 const unsigned int j) {
508 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
515 for (
unsigned int i = 0; i < n; ++i)
516 for (
unsigned int j = 0; j < m; ++j)
517 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
537 template <
typename PayloadBlockType>
553 template <
typename... Args>
315 namespace BlockLinearOperatorImplementation {
…}
579template <
typename Range,
581 typename BlockPayload,
582 typename BlockMatrixType>
590 BlockPayload(block_matrix, block_matrix)};
592 return_op.
n_block_rows = [&block_matrix]() ->
unsigned int {
593 return block_matrix.n_block_rows();
596 return_op.n_block_cols = [&block_matrix]() ->
unsigned int {
597 return block_matrix.n_block_cols();
600 return_op.block = [&block_matrix](
unsigned int i,
604 const unsigned int m = block_matrix.n_block_rows();
605 const unsigned int n = block_matrix.n_block_cols();
610 return BlockType(block_matrix.block(i, j));
613 populate_linear_operator_functions(return_op);
646template <std::size_t m,
650 typename BlockPayload>
653 const std::array<std::array<
LinearOperator<
typename Range::BlockType,
654 typename Domain::BlockType,
655 typename BlockPayload::BlockType>,
659 static_assert(m > 0 && n > 0,
660 "a blocked LinearOperator must consist of at least one block");
668 return_op.
n_block_rows = []() ->
unsigned int {
return m; };
670 return_op.n_block_cols = []() ->
unsigned int {
return n; };
672 return_op.block = [ops](
unsigned int i,
unsigned int j) ->
BlockType {
679 populate_linear_operator_functions(return_op);
700template <
typename Range = BlockVector<
double>,
701 typename Domain = Range,
702 typename BlockPayload =
703 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
704 typename BlockMatrixType>
712 BlockPayload(block_matrix, block_matrix)};
714 return_op.
n_block_rows = [&block_matrix]() ->
unsigned int {
715 return block_matrix.n_block_rows();
718 return_op.n_block_cols = [&block_matrix]() ->
unsigned int {
719 return block_matrix.n_block_cols();
722 return_op.block = [&block_matrix](
unsigned int i,
726 const unsigned int m = block_matrix.n_block_rows();
727 const unsigned int n = block_matrix.n_block_cols();
733 return BlockType(block_matrix.block(i, j));
738 populate_linear_operator_functions(return_op);
761template <std::
size_t m,
typename Range,
typename Domain,
typename BlockPayload>
765 typename Domain::BlockType,
766 typename BlockPayload::BlockType>,
770 m > 0,
"a blockdiagonal LinearOperator must consist of at least one block");
775 std::array<std::array<BlockType, m>, m> new_ops;
781 for (
unsigned int i = 0; i < m; ++i)
782 for (
unsigned int j = 0; j < m; ++j)
786 new_ops[i][j] = ops[i];
793 new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
796 return block_operator<m, m, Range, Domain>(new_ops);
810template <std::
size_t m,
typename Range,
typename Domain,
typename BlockPayload>
814 typename Domain::BlockType,
815 typename BlockPayload::BlockType> &op)
818 "a blockdiagonal LinearOperator must consist of at least "
823 std::array<BlockType, m> new_ops;
871template <
typename Range = BlockVector<
double>,
872 typename Domain = Range,
873 typename BlockPayload =
874 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
881 typename BlockPayload::BlockType(diagonal_inverse)};
901 diagonal_inverse.
block(0, 0).vmult(v.block(0), u.block(0));
902 for (
unsigned int i = 1; i < m; ++i)
904 auto &dst = v.block(i);
907 for (
unsigned int j = 0; j < i; ++j)
910 diagonal_inverse.
block(i, i).vmult(dst,
915 return_op.vmult_add = [
block_operator, diagonal_inverse](Range &v,
934 diagonal_inverse.
block(0, 0).vmult_add(v.block(0), u.block(0));
936 for (
unsigned int i = 1; i < m; ++i)
938 diagonal_inverse.
block(i, i).reinit_range_vector(
942 for (
unsigned int j = 0; j < i; ++j)
945 diagonal_inverse.
block(i, i).vmult_add(v.block(i), *tmp);
988template <
typename Range = BlockVector<
double>,
989 typename Domain = Range,
990 typename BlockPayload =
991 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
998 typename BlockPayload::BlockType(diagonal_inverse)};
1018 diagonal_inverse.
block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1020 for (
int i = m - 2; i >= 0; --i)
1022 auto &dst = v.block(i);
1025 for (
unsigned int j = i + 1; j < m; ++j)
1028 diagonal_inverse.
block(i, i).vmult(dst,
1033 return_op.vmult_add = [
block_operator, diagonal_inverse](Range &v,
1051 diagonal_inverse.
block(m - 1, m - 1)
1052 .vmult_add(v.block(m - 1), u.block(m - 1));
1054 for (
int i = m - 2; i >= 0; --i)
1056 diagonal_inverse.
block(i, i).reinit_range_vector(
1060 for (
unsigned int j = i + 1; j < m; ++j)
1063 diagonal_inverse.
block(i, i).vmult_add(v.block(i), *tmp);
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
BlockLinearOperator(const std::array< BlockType, m > &ops)
LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > BlockType
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
std::function< BlockType(unsigned int, unsigned int)> block
BlockLinearOperator(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const BlockLinearOperator< Range, Domain, BlockPayload > &)=default
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
std::function< unsigned int()> n_block_cols
std::function< unsigned int()> n_block_rows
BlockLinearOperator(const Op &op)
BlockLinearOperator(const std::array< std::array< BlockType, n >, m > &ops)
BlockLinearOperator(const BlockPayload &payload)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
EmptyBlockPayload(const Args &...)
PayloadBlockType BlockType
#define DEAL_II_NAMESPACE_OPEN
constexpr bool running_in_debug_mode()
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
static ::ExceptionBase & ExcMessage(std::string arg1)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const std::array< std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, n >, m > &ops)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &ops)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &block_operator, const BlockLinearOperator< Domain, Range, BlockPayload > &diagonal_inverse)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const BlockMatrixType &block_matrix)
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType > &op)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &block_matrix)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
void populate_linear_operator_functions(::BlockLinearOperator< Range, Domain, BlockPayload > &op)
void apply_with_intermediate_storage(const Function1 &first_op, const Function2 &loop_op, Range &v, const Domain &u, bool add)
static bool equal(const T *p1, const T *p2)