16 #ifndef dealii_block_linear_operator_h 17 #define dealii_block_linear_operator_h 19 #include <deal.II/base/config.h> 21 #include <deal.II/base/exceptions.h> 23 #include <deal.II/lac/linear_operator.h> 26 DEAL_II_NAMESPACE_OPEN
32 namespace BlockLinearOperatorImplementation
34 template <
typename PayloadBlockType =
40 template <
typename Number>
43 template <
typename Range = BlockVector<
double>,
44 typename Domain = Range,
45 typename BlockPayload =
46 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
49 template <
typename Range = BlockVector<
double>,
50 typename Domain = Range,
51 typename BlockPayload =
52 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
53 typename BlockMatrixType>
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>,
73 typename Domain = Range,
74 typename BlockPayload =
79 typename Domain::BlockType,
80 typename BlockPayload::BlockType>,
85 typename Domain = Range,
86 typename BlockPayload =
91 typename Domain::BlockType,
92 typename BlockPayload::BlockType> &op);
101 template <
typename Range = BlockVector<
double>,
102 typename Domain = Range,
103 typename BlockPayload =
104 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>,
105 typename BlockMatrixType>
109 template <
typename Range = BlockVector<
double>,
110 typename Domain = Range,
111 typename BlockPayload =
112 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
118 template <
typename Range = BlockVector<
double>,
119 typename Domain = Range,
120 typename BlockPayload =
121 internal::BlockLinearOperatorImplementation::EmptyBlockPayload<>>
200 template <
typename Range,
typename Domain,
typename BlockPayload>
202 :
public LinearOperator<Range, Domain, typename BlockPayload::BlockType>
206 typename Domain::BlockType,
207 typename BlockPayload::BlockType>;
218 typename BlockPayload::
BlockType(payload, payload))
220 n_block_rows = []() ->
unsigned int {
224 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_rows called"));
228 n_block_cols = []() ->
unsigned int {
232 "Uninitialized BlockLinearOperator<Range, Domain>::n_block_cols called"));
236 block = [](
unsigned int,
unsigned int) ->
BlockType {
240 "Uninitialized BlockLinearOperator<Range, Domain>::block called"));
256 template <
typename Op>
259 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
267 template <
size_t m,
size_t n>
270 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
281 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
294 template <
typename Op>
298 *
this = block_operator<Range, Domain, BlockPayload, Op>(op);
307 template <
size_t m,
size_t n>
309 operator=(
const std::array<std::array<BlockType, n>, m> &ops)
311 *
this = block_operator<m, n, Range, Domain, BlockPayload>(ops);
324 *
this = block_diagonal_operator<m, Range, Domain, BlockPayload>(ops);
345 std::function<BlockType(unsigned int, unsigned int)>
block;
351 namespace BlockLinearOperatorImplementation
358 template <
typename Function1,
363 apply_with_intermediate_storage(
const Function1 &first_op,
364 const Function2 &loop_op,
372 tmp->reinit(v,
true);
374 const unsigned int n = u.n_blocks();
375 const unsigned int m = v.n_blocks();
377 for (
unsigned int i = 0; i < m; ++i)
379 first_op(*tmp, u, i, 0);
380 for (
unsigned int j = 1; j < n; ++j)
381 loop_op(*tmp, u, i, j);
392 template <
typename Range,
typename Domain,
typename BlockPayload>
394 populate_linear_operator_functions(
404 for (
unsigned int i = 0; i < m; ++i)
405 op.
block(i, 0).reinit_range_vector(v.block(i), omit_zeroing_entries);
417 for (
unsigned int i = 0; i < n; ++i)
418 op.
block(0, i).reinit_domain_vector(v.block(i), omit_zeroing_entries);
423 op.
vmult = [&op](Range &v,
const Domain &u) {
431 const auto first_op = [&op](Range & v,
433 const unsigned int i,
434 const unsigned int j) {
435 op.
block(i, j).vmult(v.block(i), u.block(j));
438 const auto loop_op = [&op](Range & v,
440 const unsigned int i,
441 const unsigned int j) {
442 op.
block(i, j).vmult_add(v.block(i), u.block(j));
445 apply_with_intermediate_storage(first_op, loop_op, v, u,
false);
449 for (
unsigned int i = 0; i < m; ++i)
451 op.
block(i, 0).vmult(v.block(i), u.block(0));
452 for (
unsigned int j = 1; j < n; ++j)
453 op.
block(i, j).vmult_add(v.block(i), u.block(j));
458 op.
vmult_add = [&op](Range &v,
const Domain &u) {
466 const auto first_op = [&op](Range & v,
468 const unsigned int i,
469 const unsigned int j) {
470 op.
block(i, j).vmult(v.block(i), u.block(j));
473 const auto loop_op = [&op](Range & v,
475 const unsigned int i,
476 const unsigned int j) {
477 op.
block(i, j).vmult_add(v.block(i), u.block(j));
480 apply_with_intermediate_storage(first_op, loop_op, v, u,
true);
484 for (
unsigned int i = 0; i < m; ++i)
485 for (
unsigned int j = 0; j < n; ++j)
486 op.
block(i, j).vmult_add(v.block(i), u.block(j));
490 op.
Tvmult = [&op](Domain &v,
const Range &u) {
498 const auto first_op = [&op](Range & v,
500 const unsigned int i,
501 const unsigned int j) {
502 op.
block(j, i).Tvmult(v.block(i), u.block(j));
505 const auto loop_op = [&op](Range & v,
507 const unsigned int i,
508 const unsigned int j) {
509 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
512 apply_with_intermediate_storage(first_op, loop_op, v, u,
false);
516 for (
unsigned int i = 0; i < n; ++i)
518 op.
block(0, i).Tvmult(v.block(i), u.block(0));
519 for (
unsigned int j = 1; j < m; ++j)
520 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
525 op.
Tvmult_add = [&op](Domain &v,
const Range &u) {
533 const auto first_op = [&op](Range & v,
535 const unsigned int i,
536 const unsigned int j) {
537 op.
block(j, i).Tvmult(v.block(i), u.block(j));
540 const auto loop_op = [&op](Range & v,
542 const unsigned int i,
543 const unsigned int j) {
544 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
547 apply_with_intermediate_storage(first_op, loop_op, v, u,
true);
551 for (
unsigned int i = 0; i < n; ++i)
552 for (
unsigned int j = 0; j < m; ++j)
553 op.
block(j, i).Tvmult_add(v.block(i), u.block(j));
574 template <
typename PayloadBlockType>
575 class EmptyBlockPayload
590 template <
typename... Args>
616 template <
typename Range,
618 typename BlockPayload,
619 typename BlockMatrixType>
627 BlockPayload(block_matrix, block_matrix));
629 return_op.
n_block_rows = [&block_matrix]() ->
unsigned int {
630 return block_matrix.n_block_rows();
633 return_op.
n_block_cols = [&block_matrix]() ->
unsigned int {
634 return block_matrix.n_block_cols();
637 return_op.
block = [&block_matrix](
unsigned int i,
638 unsigned int j) -> BlockType {
640 const unsigned int m = block_matrix.n_block_rows();
641 const unsigned int n = block_matrix.n_block_cols();
646 return BlockType(block_matrix.block(i, j));
649 populate_linear_operator_functions(return_op);
686 typename BlockPayload>
689 const std::array<std::array<
LinearOperator<
typename Range::BlockType,
690 typename Domain::BlockType,
691 typename BlockPayload::BlockType>,
695 static_assert(m > 0 && n > 0,
696 "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 {
716 populate_linear_operator_functions(return_op);
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));
773 populate_linear_operator_functions(return_op);
796 template <
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 <
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)
941 block_operator.
block(i, j).vmult_add(dst, v.block(j));
943 diagonal_inverse.
block(i, i).vmult(dst,
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)
976 block_operator.
block(i, j).vmult_add(*tmp, v.block(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)
1057 block_operator.
block(i, j).vmult_add(dst, v.block(j));
1059 diagonal_inverse.
block(i, i).vmult(dst,
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)
1092 block_operator.
block(i, j).vmult_add(*tmp, v.block(j));
1094 diagonal_inverse.
block(i, i).vmult_add(v.block(i), *tmp);
1103 DEAL_II_NAMESPACE_CLOSE
BlockLinearOperator(const std::array< BlockType, m > &ops)
BlockLinearOperator(const Op &op)
std::function< void(Range &v, bool omit_zeroing_entries)> reinit_range_vector
std::function< void(Domain &v, const Range &u)> Tvmult_add
std::function< void(Domain &v, const Range &u)> Tvmult
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< BlockType, m > &ops)
BlockLinearOperator< Range, Domain, BlockPayload > block_operator(const BlockMatrixType &matrix)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const Op &op)
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
LinearOperator< Range, Domain, Payload > null_operator(const LinearOperator< Range, Domain, Payload > &)
BlockLinearOperator< Range, Domain, BlockPayload > & operator=(const std::array< std::array< BlockType, n >, m > &ops)
static::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
static::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
std::function< void(Range &v, const Domain &u)> vmult
std::function< void(Domain &v, bool omit_zeroing_entries)> reinit_domain_vector
BlockLinearOperator< Range, Domain, BlockPayload > block_diagonal_operator(const std::array< LinearOperator< typename Range::BlockType, typename Domain::BlockType, typename BlockPayload::BlockType >, m > &)
BlockLinearOperator(const std::array< std::array< BlockType, n >, m > &ops)
std::function< void(Range &v, const Domain &u)> vmult_add
std::function< BlockType(unsigned int, unsigned int)> block
std::function< unsigned int()> n_block_cols
BlockLinearOperator(const BlockPayload &payload)
EmptyBlockPayload(const Args &...)
std::function< unsigned int()> n_block_rows
PayloadBlockType BlockType
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_back_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)
static bool equal(const T *p1, const T *p2)
LinearOperator< Domain, Range, typename BlockPayload::BlockType > block_forward_substitution(const BlockLinearOperator< Range, Domain, BlockPayload > &, const BlockLinearOperator< Domain, Range, BlockPayload > &)