16 #ifndef dealii_block_linear_operator_h
17 #define dealii_block_linear_operator_h
32 namespace BlockLinearOperatorImplementation
34 template <
typename PayloadBlockType =
36 class EmptyBlockPayload;
40 template <
typename Number>
43 template <
typename Range = BlockVector<
double>,
44 typename Domain = Range,
45 typename BlockPayload =
46 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
50 template <
typename Range = BlockVector<
double>,
51 typename Domain = Range,
52 typename BlockPayload =
53 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
54 typename BlockMatrixType>
58 template <std::size_t m,
61 typename Domain = Range,
62 typename BlockPayload =
66 const std::array<std::array<
LinearOperator<
typename Range::BlockType,
67 typename Domain::BlockType,
68 typename BlockPayload::BlockType>,
72 template <std::size_t m,
74 typename Domain = Range,
75 typename BlockPayload =
80 typename Domain::BlockType,
81 typename BlockPayload::BlockType>,
84 template <std::size_t m,
86 typename Domain = Range,
87 typename BlockPayload =
92 typename Domain::BlockType,
93 typename BlockPayload::BlockType> &op);
102 template <
typename Range = BlockVector<
double>,
103 typename Domain = Range,
104 typename BlockPayload =
105 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
106 typename BlockMatrixType>
110 template <
typename Range = BlockVector<
double>,
111 typename Domain = Range,
112 typename BlockPayload =
113 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
119 template <
typename Range = BlockVector<
double>,
120 typename Domain = Range,
121 typename BlockPayload =
122 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
201 template <
typename Range,
typename Domain,
typename BlockPayload>
203 :
public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
207 typename Domain::BlockType,
208 typename BlockPayload::BlockType>;
219 typename BlockPayload::
BlockType(payload, payload))
225 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
233 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
241 "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
257 template <
typename Op>
260 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
268 template <std::
size_t m, std::
size_t n>
271 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
279 template <std::
size_t m>
282 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
295 template <
typename Op>
299 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
308 template <std::
size_t m, std::
size_t n>
310 operator=(
const std::array<std::array<BlockType, n>, m> &ops)
312 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
321 template <std::
size_t m>
325 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
352 namespace BlockLinearOperatorImplementation
359 template <
typename Function1,
365 const Function2 &loop_op,
373 tmp->reinit(v,
true);
375 const unsigned int n = u.n_blocks();
376 const unsigned int m = v.n_blocks();
378 for (
unsigned int i = 0; i < m; ++i)
380 first_op(*tmp, u, i, 0);
381 for (
unsigned int j = 1; j < n; ++j)
382 loop_op(*tmp, u, i, j);
393 template <
typename Range,
typename Domain,
typename BlockPayload>
405 for (
unsigned int i = 0; i < m; ++i)
406 op.
block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
418 for (
unsigned int i = 0; i < n; ++i)
419 op.
block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
424 op.
vmult = [&op](Range &v,
const Domain &u) {
432 const auto first_op = [&op](Range & v,
434 const unsigned int i,
435 const unsigned int j) {
436 op.
block(i, j).vmult(v.block(i), u.block(j));
439 const auto loop_op = [&op](Range & v,
441 const unsigned int i,
442 const unsigned int j) {
443 op.
block(i, j).vmult_add(v.block(i), u.block(j));
450 for (
unsigned int i = 0; i < m; ++i)
452 op.
block(i, 0).vmult(v.block(i), u.block(0));
453 for (
unsigned int j = 1; j < n; ++j)
454 op.
block(i, j).vmult_add(v.block(i), u.block(j));
459 op.
vmult_add = [&op](Range &v,
const Domain &u) {
467 const auto first_op = [&op](Range & v,
469 const unsigned int i,
470 const unsigned int j) {
471 op.
block(i, j).vmult(v.block(i), u.block(j));
474 const auto loop_op = [&op](Range & v,
476 const unsigned int i,
477 const unsigned int j) {
478 op.
block(i, j).vmult_add(v.block(i), u.block(j));
485 for (
unsigned int i = 0; i < m; ++i)
486 for (
unsigned int j = 0; j < n; ++j)
487 op.
block(i, j).vmult_add(v.block(i), u.block(j));
491 op.
Tvmult = [&op](Domain &v,
const Range &u) {
499 const auto first_op = [&op](Range & v,
501 const unsigned int i,
502 const unsigned int j) {
503 op.
block(j, i).Tvmult(v.block(i), u.block(j));
506 const auto loop_op = [&op](Range & v,
508 const unsigned int i,
509 const unsigned int j) {
510 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
517 for (
unsigned int i = 0; i < n; ++i)
519 op.
block(0, i).Tvmult(v.block(i), u.block(0));
520 for (
unsigned int j = 1; j < m; ++j)
521 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
526 op.
Tvmult_add = [&op](Domain &v,
const Range &u) {
534 const auto first_op = [&op](Range & v,
536 const unsigned int i,
537 const unsigned int j) {
538 op.
block(j, i).Tvmult(v.block(i), u.block(j));
541 const auto loop_op = [&op](Range & v,
543 const unsigned int i,
544 const unsigned int j) {
545 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
552 for (
unsigned int i = 0; i < n; ++i)
553 for (
unsigned int j = 0; j < m; ++j)
554 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
575 template <
typename PayloadBlockType>
591 template <
typename... Args>
617 template <
typename Range,
619 typename BlockPayload,
620 typename BlockMatrixType>
628 BlockPayload(block_matrix, block_matrix)};
630 return_op.
n_block_rows = [&block_matrix]() ->
unsigned int {
631 return block_matrix.n_block_rows();
634 return_op.n_block_cols = [&block_matrix]() ->
unsigned int {
635 return block_matrix.n_block_cols();
638 return_op.block = [&block_matrix](
unsigned int i,
639 unsigned int j) -> BlockType {
641 const unsigned int m = block_matrix.n_block_rows();
642 const unsigned int n = block_matrix.n_block_cols();
647 return BlockType(block_matrix.block(i, j));
683 template <std::size_t m,
687 typename BlockPayload>
690 const std::array<std::array<
LinearOperator<
typename Range::BlockType,
691 typename Domain::BlockType,
692 typename BlockPayload::BlockType>,
696 static_assert(m > 0 && n > 0,
697 "a blocked LinearOperator must consist of at least one block");
705 return_op.
n_block_rows = []() ->
unsigned int {
return m; };
707 return_op.n_block_cols = []() ->
unsigned int {
return n; };
709 return_op.block = [ops](
unsigned int i,
unsigned int j) -> BlockType {
737 template <
typename Range,
739 typename BlockPayload,
740 typename BlockMatrixType>
748 BlockPayload(block_matrix, block_matrix)};
750 return_op.
n_block_rows = [&block_matrix]() ->
unsigned int {
751 return block_matrix.n_block_rows();
754 return_op.n_block_cols = [&block_matrix]() ->
unsigned int {
755 return block_matrix.n_block_cols();
758 return_op.block = [&block_matrix](
unsigned int i,
759 unsigned int j) -> BlockType {
761 const unsigned int m = block_matrix.n_block_rows();
762 const unsigned int n = block_matrix.n_block_cols();
768 return BlockType(block_matrix.block(i, j));
796 template <std::
size_t m,
typename Range,
typename Domain,
typename BlockPayload>
800 typename Domain::BlockType,
801 typename BlockPayload::BlockType>,
805 m > 0,
"a blockdiagonal LinearOperator must consist of at least one block");
810 std::array<std::array<BlockType, m>, m> new_ops;
816 for (
unsigned int i = 0; i < m; ++i)
817 for (
unsigned int j = 0; j < m; ++j)
821 new_ops[i][j] = ops[i];
828 new_ops[i][j].reinit_domain_vector = ops[j].reinit_domain_vector;
831 return block_operator<m, m, Range, Domain>(new_ops);
845 template <std::
size_t m,
typename Range,
typename Domain,
typename BlockPayload>
849 typename Domain::BlockType,
850 typename BlockPayload::BlockType> &op)
853 "a blockdiagonal LinearOperator must consist of at least "
858 std::array<BlockType, m> new_ops;
907 template <
typename Range,
typename Domain,
typename BlockPayload>
914 typename BlockPayload::BlockType(diagonal_inverse)};
934 diagonal_inverse.
block(0, 0).vmult(v.block(0), u.block(0));
935 for (
unsigned int i = 1; i < m; ++i)
937 auto &dst = v.block(i);
940 for (
unsigned int j = 0; j < i; ++j)
943 diagonal_inverse.
block(i, i).vmult(dst,
948 return_op.vmult_add = [
block_operator, diagonal_inverse](Range & v,
967 diagonal_inverse.
block(0, 0).vmult_add(v.block(0), u.block(0));
969 for (
unsigned int i = 1; i < m; ++i)
971 diagonal_inverse.
block(i, i).reinit_range_vector(
975 for (
unsigned int j = 0; j < i; ++j)
978 diagonal_inverse.
block(i, i).vmult_add(v.block(i), *tmp);
1022 template <
typename Range,
typename Domain,
typename BlockPayload>
1029 typename BlockPayload::BlockType(diagonal_inverse)};
1049 diagonal_inverse.
block(m - 1, m - 1).vmult(v.block(m - 1), u.block(m - 1));
1051 for (
int i = m - 2; i >= 0; --i)
1053 auto &dst = v.block(i);
1056 for (
unsigned int j = i + 1; j < m; ++j)
1059 diagonal_inverse.
block(i, i).vmult(dst,
1064 return_op.vmult_add = [
block_operator, diagonal_inverse](Range & v,
1082 diagonal_inverse.
block(m - 1, m - 1)
1083 .vmult_add(v.block(m - 1), u.block(m - 1));
1085 for (
int i = m - 2; i >= 0; --i)
1087 diagonal_inverse.
block(i, i).reinit_range_vector(
1091 for (
unsigned int j = i + 1; j < m; ++j)
1094 diagonal_inverse.
block(i, i).vmult_add(v.block(i), *tmp);